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A Numerical Study of the Relative Class 
Numbers of Real Quadratic Integral Domains 


By Harvey Cohn 


1. Introduction. In a classic paper in 1856 Dirichlet gave some applications of a 
formula for the ratio of the class number of a quadratic integral domain in a real 
field to the class number of the whole integral domain (of all quadratic integers in 
that field), with the principal objective of showing that this ratio takes many 
values (such as 1) infinitely often for the real case, in support of a conjecture of 
Gauss. 

The object of this paper is first of all to give Dirichlet’s results briefly, together 
‘with some theorems and illustrations immediately deducible from them (in order 
to restrict the computation to cases in which the theory is of more help). We shall, 
of course, offer various tables of relative class numbers, such data being our main 
object. We emphasize quadratic integral domains of prime power conductor under 
the whole integral domain (of all quadratic integers of the field). 

Weask, in particular, when the relative class number is divisible by 2 and 4, and 
find simple linear congruence conditions. When we ask which prime conductors 
have relative class numbers divisible by 3, we find such primes are essentially the 
splitting primes of certain cubic fields and therefore representable by quadratic 
forms, according to the classic work of Dedekind [3]. This is basically an application 

| of class-field theory and perhaps the tables emerging would be of some experimental 
use. The classic background is amplified in [7], [5], and [2]. 

Here it might be appropriate to remark that the tables-given below have a 
“natural” limit of diminishing returns owing to the fact that the relevant portions 
of classical algebraic number theory were developed long ago with relatively little 
data, and it would be desirable to see the theory profit from more data before 
great feats of computer endurance are attempted. 


2. Notation and Terminology. We follow the convention that Latin letters 
generally denote rational integers and Greek letters denote algebraic integers. The 
following symbols and terms appear throughout the work: 


m is a square-free integer > 1. 
R(m*”) is the field generated by m*”. 
is the discriminant of the field generated by m'?.d = mifm = 1 
(mod 4), d = 4m if m # 1 (mod 4). 
is the indicator of the type of field, c = 2 if m = 1 (mod 4),c = 1 
if m = 1 (mod 4). Thus d = 4m/c’. 
is the set of all algebraic integers of R(m™*). It consists of w = 
(x + ym"*)/c for which z and y are rational integers subject 
only to the condition x = y (mod c). 


Received August 18, 1961. This research was supported in part by a National Science 
Foundation grant. Computer service on GEORGE was donated by the Argonne National 
Laboratories of the U. S. Atomic Energy Commission through the cooperation of Dr. William 

F. Miller, Director of the Applied Mathematics Division. 
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F(w) is the function defined by y = F(w) in the above definition. 

Dy is an arbitrary integral domain (ring with unity) in R(m””), given 
uniquely for any integer f > 0. It consists of the subset of alge- 
braic integers w in © for which f| F(w). Here f is called the 
conductor. It is the index of D; in D, and OD; = ©. 

frd is the ring-discriminant of D; . Its purpose is that any given D(>1) 
which is = 0 or 1 (mod 4) can be written uniquely as f'd for some 
f and d. Thus f‘d completely determines R(m™”) and D; in R(m*”). 


h(f'd) is the class number (of ideals prime to f) in D; . 

H(f) is the relative class number = h(f'd)/h(d), (used when the value 
of d, or the field, is understood in context). 

€ is the fundamental unit, written « = (a + bm””)/c. Here a > 0, 
b > O and a = b (mod ce). 

e is the norm of the fundamental unit, actually +1, N(e) = e = 
(a? — mb’)/c’. 

¥v(f) is the value of f11(1 — (d/q)/q) extended over primes q which divide 
f. Here (d/q) is the Kronecker residue symbol, (thus (d/2) = 
(2/d)). 

o(f) is the minimum exponent ¢ (>1) for which ¢ ¢ OD; or for which 


f|F(é). It can be shown directly that ¢(f) | ¥(f). By classical 
methods of primitive root theory, if f| F(e*), then ¢(f) | u. 


Other symbols appear only locally and can best be defined as they arise. 


3. Dirichlet’s Theorems. The starting point is the following theorem, in prin- 
ciple due to Gauss: For a given field R(m*™*), (with m > 0), 


(3.1) H(f) = V(f)/o(h). 


If m < 0, the formula is modified so that, for instance, with f > 1, ¢(f) is replaced 
by half the number of units in ©. (We do not need the modified formula for the 
machine part of the calculation, but for supporting computations in Section 7). 

Now y(f) is fairly easy to find, but the calculation of ¢(f) is the part requiring 
the electronic computer. Dirichlet [4] showed, however, that if f = p,"" --- p,”*, 
where the primes p; come from a given finite set, then the values of H(f) also come 
from a finite set as the exponents F; vary; in fact H(f) = Hy, a constant if each 
F; is sufficiently large. An examination of Dirichlet’s method leads to the rule that 
if p; is odd and fy is such that H(fo) = H(fop:) (while if one p; = 2, fo satisfies 
A (fo) = H(4fo)), then H(f) = H(fo) if fo| f, (recalling the prime divisors of f 
are to be limited to the p,). 

From general principles it also follows that if f |g, then H(f) | H(g). 

The main step in understanding these results is to consider any f which contains 
all the odd primes p; (and possibly 2”) as divisors. Then f | F(e*”), ie. &” = 
(xs + yym"™) /c, where f | yy . But, let f* be the factor of y, consisting of powers of 
the p;. (Thus f | f* while (y,/f*, f) = 1.) Then for p; odd, F(e&””') = p.f*g, 
where (g, f) = 1, as we prove by using the binomial theorem, (in a manner reminis- 
cent of the proof that a primitive root modulo p’ is a primitive root modulo p”, 
n > 2). For p; even, special attention must be given the denominator c = 2, but 
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this can be left to the reader, as well as the completion of the proof of the above 
results by induction. 

If we restrict f to powers of a prime p then we find H(p"**) = H(p”") or pH(p") 
(n = 1), but eventually H(p"™) = H(p") then H(p”) = H(p") for allm = n 
when p is odd, while for p = 2, H(2"**) = H(2") or 2H(2"), (n = 1), but even- 
tually H(2"**) = H(2") where H(2”) = H(2") for all n = m. 

4. Simple Cases. We first consider those g which divide 6 mb. In these cases the 
values of H(q‘) are easily seen by elementary hand calculations, and we often omit 
these from the tables to make room for more interesting values. 


f=2" 
(1 if = 1, 
Define M(a, b,f) =4min (2",f) if 24 =1, 


f 
\min (2*,f/2) if 2*>1, f 


IV 


9 


‘ 
aH, 


IV 


where 2% || a, i.e., 2* | a but 2“** 4 a, and likewise 2” || b. Then if d = 0 (mod 4), 
(4.1) H(f) = M(a,b,f), 
while if d = 1 (mod 4) and 2| ab, 
(4.2) H(f) = [2 + (d/2)|M(a/2, b/2, f/2), 
and if d = 1 (mod 4) and 2 ~ ab (whence d = 5 (mod 8)), 
(4.3) H(f) = M({a° — 3¢)/2, [a° — e]/2, f/2). 
(Note that ({a + bm’*|/2)* = ala’ — 3e]/2 + bla’ — ejm*?/2). 
f = 3” 


Let 3° Il b, 3“ || a. If 3| m, let 3° || 3a” + mb’, then 
(2 min (f,3°) if 3% =1 
(4.4) H(f) = 4 . a hen 
min (f,3°) if 3 > 1. 
If, however, 3 ~ m, let 37 || a + b’m, then 
(2 min (f, 37) if 3+ ab, 
(4.5) H(f) = <43{1 — (d/3)/3] min (f,3*) if 3a, 
\{1 — (d/3)/3] min (f,3”) if 3|6. 


~ oF : . ; . . 5 oF ‘ Fr 
(Note that f = 3” is “special” because of consideration of 3°. Compare f = q 
below). 


f=q" 
Here let g be a prime ~2, 3 for which q | mb, and let q” || b. Then 
(4.6) H(f) = min (f, 9”). 
Thus in many cases where q | m and q + 6b, then H(q") = 1 forall n, giving the 
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easiest illustration of Dirichlet’s original objective; e.g., for m = d 
for all n > 0. 


Il 
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or 
3 

— 

ll 
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5. The Program. The basic sub-routine considers the input 
(5.1) m, a,b, f 


from which ¢(f), ¥(f), and H(f) are calculated. The machine forms by induction 
e' = [a(t) + b(t)m'’|/c stored as a(t), b(t) calculated modulo f*. Then letting 
t = 1, 2, the machine records the earliest t{/= ¢(f)] for which b(t) = 0 (mod f). 
The machine next calculates ¥(f) by examining the prime factors q of f sequen- 
tially. The machine finds (d/q) for g + 2d by actually testing the solvability in 
x of x” = d (mod q), while for q | 2d, (d/q) is determined directly from the rules. 
Finally, H(f) = W(f)/o(f). The output for each input consists of 


(5.2) f, H(f), Wf), b@(f))/f (mod f). 


The last value is desired for purposes of testing F(<°”’). For example, if 


((6(f))/F,f) = 1, 
then f is a suitable fy for Section 3. 

The basic sub-routine was used in several ways. 

In one run the basic sub-routine was set up to increment f by 1 automatically 
over a range fi S f S fe where f; and f2 are given in addition to the initial data. 
For m = 5 the problem was run up to f = 4400 and for m = 2 and 3, it was run up 
to f = 1000. 

In another variation, the values of f were incremented as before but were re- 
stricted to primes in the preassigned range. (We always use the letter p to denote 
a prime.) These main runs were made for f = p an odd prime up to 997 for 38 
values of m, namely 


(5.2) Series A: 2 


IIA 


square free m = 1 (mod 4) 


IIA 


42 
(5.3) Series B: 5 


lA 


square free m = 1 (mod 4) S 97. 


The problem was programmed for the GEORGE computer with only approxi- 
mately 500 words of a 4096-word high-speed memory involved. The machine is 
internally binary with 40-bit word length and approximate speed of 50,000 two- 
address operations per second. 

In all the runs, the output consisted of the input data (5.1) (as a heading) fol- 
lowed by the output data (5.2) listed “on-line” (parallel) with the computation. 
The input and output were in decimal (internally converted) and on paper tape 
originally (but the output was later transformed to magnetic tape just to speed up 
the printing process from flexowriter to line printer). The actual input and output 
times were negligible. 

The running time for each case was about f/50 seconds. The calculations were 
run between December 1960 and May 1961. 


6. Use of Some Cyclic Groups. Let m be given and let p + 2m be an arbitrary 
given prime. Define a group in which the elements a, are the following sets: 


(6.1) a, = {x2 + ym}, where x =ty and N(x + ym’) = 0 (mod Pp), 
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and 
(6.2) a,j = {x}, where z <0. 


The group operation is multiplication (mod p), easily shown to be independent of 
the representative. When (m/p) = —1, there are p + 1 of these elements, while 
when (m/p) = +1 there are p — 1 of these elements (by excluding two values of 
t for which f = m (mod p)). In general, we have a group &, with Pp - 
(m/p) = ¥(p) elements, and with a,, as the unit element. 

We see that the group 4, is cyclic. This is true where (m/p) = —1 since the 
group is a sub-group of the cyclic group of reduced residues of algebraic integers 
modulo p, (now an ideal prime). When (m/p) = —1 we rewrite a, = alu] where 


(6.3) alu] = {x(r(1 + u) + m’ *(1 — u))}. 


Here r satisfies r° = m (mod p) and ¢ and uw are related by t = r(i + u)/(1 — w) 
(mod p). We can verify afujalv] = a[uv], hence when (m/p) = 1, Y, is isomorphic 
to the multiplicative (cyclic) residue group of rational integers modulo p. 

The important result for us is the following: if p + 2m and if r is a given integer 
dividing p — (m/p) a necessary and sufficient condition that r| H(p) is that 
ce belong to an a, which is an r-th power in %, . This result follows from the cyclic 
structure of %, once we note that (ce)*”’ = z (mod p) for z an integer, hence 
(ce)*” belongs to a, the unit element, while ¥(p) is the order of the group. 

For illustration, we start with r = 2, and take p + 2mb. Set 


a+ bm? = (x + ym *\k, or, 


(a = k(a” + ym) 
(6.4) 
\b = 2kry. 
This system is solvable, for k ¥ 0, if and only if the equation 
(6.5) ba” — 2ary + bmy’ = 0 mod p 
is solvable, with (x, y) + (0,0). The discriminant is 4c’e. Hence if N(«) = e = —1, 
then 2 | H(p), for p ~ 2mb. 
Thus for some cases, e.g., where N(e) = +1, the only possible f for which 
H(f) = 1 must come from primes in the special cases in Section 4 above. (We 


recall that if f | g, then H(f) | H(g)). Thus for m = 3, the only f for which H(f) = 1 
are now seen to be f = 3‘ and f = 2-3". 

We next consider the sub-group of 4, , called B,, all of whose elements have 
norms which are quadratic residues of p. Thus a, is necessarily in 8, , while a, 
is in B, if and only if ({f — m|/p) = +1. It is easily seen that the norms of repre- 
sentatives in a, are not all residues, by results on successions of residues and non- 
residues. Thus %, has only order (p — (m/p))/2, since it must then be of index 2. 
Now if we normalize the representative of a; in (6.1) belonging to 8, to be plus 
or minus an element of norm 1, we can say that if e = 1, then ¢ represents a perfect 
square in %, if and only if for some integers x and y 
(6.6) +c'e = (x + ym’)? (mod p). 

But the condition for a perfect square in %, is precisely the condition that +e 


represents a perfect fourth power in &%,, or 4| H(p). Expanding (6.6), we dis- 
cover we must be able to solve simultaneously 
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(tea =2° +m’ 
(6.7) | (mod p). 
|eb = 2xy 


An elementary calculation reveals this system is solvable if and only if (with signs 


8, %& = +1), 
( + my” 


2 2 
x — my = & 


8\Ca 


(6.8) (mod p). 


I 


For this it is necessary and sufficient that 2c(s;a + s.) and 2mec(s;a — 8) be per- 
fect squares modulo p. With some manipulation, we find, if N(e) = e = +1 and 
p + 2mb, then a necessary and sufficient condition that 4| H(p) is that 


(6.9) (—1/p) = (m/p) = ([2a/c — 2]/p). 

We can often simplify the result (6.9) to take the form 

(6.10) (—S/p) = (Q/p) = (R/p), 

for smaller values of Q and R shown in the columns 10 and 11 of Table I with S = 1, 
except for m = 15 and 35, where S = 2. When e = —1, there are still many oc- 


currences of H(p) = 4 (the smallest such value is listed in column 11). 


7. Divisibility by 3. A more interesting case is r = 3. This can occur (for p + 6mb) 
only when 3 | ¥(p) or (—3m/p) = 1. We ask, when can we solve ce = 
k(x + ym™*)® (mod p) or 


(a = kix* + 32y'm] 
(7.1) 4 4 ; (mod p), 

\b = kBa'y + ym) 
for zy + 0? Eliminating k, we see this leads to the solvability of A(2/y) 
= 0 mod p where d is a polynomial defining a root of a cubic field, 


(7.2) A(E) = be — Bak? + 3bEm — am = O. 
Hence 3 | H(p) (for p ¢ 6m) if and only if p is a splitting prime for the field R(£). 
In fact, p must split into three distinct prime ideals since (—3m/p) = 1, and the 


discriminant D,; of the cubic can be shown to differ from —3m by a rational square. 
The reader is referred to Hasse’s work [6] for details on the method. 

Finding the field discriminant of R(é) is rather lengthy but since the methods 
are so well-known we can merely outline the steps. The module [1, bé, am/é] con- 
sists only of integers of R(£) and its discriminant is —108mc* by a direct calcula- 
tion. Since only perfect squares could be superfluous factors of the discriminant, 
we need examine the basis elements to see if r + sbt + tam/é can be divisible by 
2 (or 3) without 7, s, and ¢ being simultaneously divisible by 2 (or 3). We find 
the only possibilities are the following cases which we leave for the reader to verify: 


Case i. 3)| mand 3 |b; then 3 | bé and 3 |(am/é) 
Case ii. 34 mand 9| a (or b); then 3 | bE (or 3 | (am/é)) 
Case iii. 3 7 mab and am = +b (mod 9); then 3 | (bE + esexam/E — e2) 


where e; = +1 = am, & = +1 = b (mod 3) 
Case iv. c = 2; then 2 | bé, 4 |(bE + am/é). 
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These calculations were made partly on the basis of possible ideal factorizations 
of (2) and (3) and partly as a direct consequence of the following equation for 
u = (bE + am/é)e: 


uw — 3b(1 — m)p’ + 3(B’ — a’) (1 — m)p 
(7.3) + [a’(6m + 2) + a’b(m’® — 12m + 3) 
— ab’(2m + 6m’) + b°(9m? — 1)] = 0. 


The occurrences of cases (i-iii) are noted in column 7 of Table I. 
We finally obtain 


(7.4) dxf; = D; = —108m/s‘c’, 

where 
\ =9 if 3|m,3\b, 

(7.5) \8 =3 if 34m,9|ab, orif 34 mab, am = +b (mod 9), 
\s = 1 otherwise. 


We then consider the set of h(d;f;’) primitive reduced quadratic forms of dis- 
criminant D;. Those which are perfect cubes under composition represent pre- 
cisely all primes p(+ 6m) for which 3 | H(p). 

A supporting computation was made by Mr. Roy Lippmann on an IBM 650 
to calculate all primitive reduced forms from D;. The square-free kernel m, is 
shown in Table I, together with h(D;) and the conductor f, . The h(D;) primitive 
forms (A, B, C) which are cubes under composition were most easily identified 
by finding some “convenient” small prime (p 7 6m) represented by the form and 
checking H(p), (see [1]). The coefficients A and B of forms and representative 
primes p and H(p) are listed in Table III. 

Now in every case, it so happens that 3 || h(d,f;’), hence there are h(dsf;) /3 
forms which are perfect cubes. Also, the ambiguous forms are always perfect cubes, 
but in general they are not the complete set. The non-ambiguous forms, naturally, 
are written two at a time by means of +B. 


8. Irregular Primes. We finally note that there are many odd primes p, for 
which, for some fixed 7 > 0, 


(8.1) H(p") = H(p) min (p™”, p’). 


We call these primes irregular and we call i the index of irregularity. When p + 6mb 
such cases are explained by some combinational curiosities much less transparent 
than those occurring in Section 3. They are listed because the occurrence of prime 
divisors of f in the relative class number is of some theoretical value. 

These values were found by scanning the outputs (5.2) as f ran over the odd 
primes p for cases where b = 0 (mod p). The 53 individual cases which emerged 
were tested by rerunning these cases, using f = p. The values of b/f + 0 (mod f) 
indicated primes of index 1, while those where b/f = 0 while b/fp = 0 (mod f) 
indicated primes p of index 2. No odd primes of higher index emerged from the 
experiment. 
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9. Summary of Calculations. The problem ran some 40 hours and generated 
some 300 pages of table obviously too much to reproduce! We therefore attempt 
a qualitative résume. 

From the output, we would readily believe that when e = —1 there are in- 
finitely many odd primes for which H(p) = 1, while when e = 1 there are in- 
finely many primes for which H(p} = 2. Indeed, even in the case e = 1, we know 
(from Section 4) that if p | m and p + 6b then H(p) = 1. In either case, except for 
scattered irregular primes in Table IV, H(p) = H(p”). 


A frequency count is surprising in its uniformity. When e = —1, we examine 
the 167 odd primes < 1000 and find H(p) = 1 in 39-43 per cent of these primes 
as m varies, while when e = +1 the corresponding case H(p) = 2 occurs for 


56-63 per cent of these primes as m varies. If we define P(m, n; x) as the proportion 
of primes S x for which H(p) = n (in reference to R(m”)) we find a reasonably 
steady value for P(5, 1; x). For instance, P(5, 1; 500) = 42 per cent, P(5, 1; 
1000) = 41 per cent, P(5, 1; 2000) = 39 per cent, P(5, 1; 4000) = 37 per cent. 

Continuing with m = 5, H(p) (as far as we might imagine) “should’’ take all 
prime values but it seems to take large values “rather slowly.”” The earliest p for 
some larger primes are H(911) = 13, H(1087) = 17, H(3079) = 19, H(1103) = 
23. For p < 4400, H(p) takes no larger prime! Thus an “asymptotic” study of 
the values of H(p) can be expected to be astronomical in size (perhaps larger 
than for studies of classical prime number distributions) . 

Table II is given to point out some relative class numbers which are small 
prime powers; H(p) = 3 is in Table III; and H(p) = 2 or 4 comes from columns 
10 and 11 of Table I. Despite the uniformity of the earlier frequency count, some 
values of m seem to be more “amenable” to given values of H(p) than others. 
This seeming paradox might again be a manifestation of the fact that ‘“p < 1000” 
is a miniscule range of values! 

As far as congruence properties of H(p) are concerned, Sections 6 and 7 provide 
us with much more guidance. For example, by the uniform density of primes in 
linear congruence classes for a fixed modulus, when e = —1, H(p) = 0 (mod 4) 
only one-third as often as H(p) = 0 (mod 4). 

In a similar manner, using known results on the distribution of primes repre- 
sented by quadratic forms [8], we can see that if k; of the h(D3) forms are perfect 
cubes, then k;/2h(D;) is the proportion of primes for which H(p) = 0 (mod 3), 
at least by ‘Dirichlet density.” Actual frequency counts show the proportion to 
be reassuringly close to 14; (with kz; = h(D3)/3 in the cases treated here). 

The congruence properties H(p) = 0 (mod 5), however, provide too few in- 
stances in the range p < 1000, to make a frequency count meaningful. 

The conditions on p which make H(p) = 0 (mod 4) when e = —1, are more 
provocative. The percentage of such p( <1000) varies from 4 per cent (when 
m = 37) to 12 per cent (when m = 89). There seems to be no simple explanation 
(e.g., in terms of linear or quadratic forms). As a matter of curiosity, when m = 5, 
H(p) = 0 (mod 4) for 


p = 61, 89, 109, 149, 269, 389, 401, 521, 661, 701, 761, 769, 809, 821, 829; 
when m = 37, this holds for 


p = 53, 101, 181, 293, 349, 397, 593; 
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TaBLeE I 
Summary of Calculation : 
Columns 1-5 are explained in Section 2, Columns 6-9 in Section 7, 
Columns 10-11 in Section 6. 























| j 
1 | 2 3 ha.S t 2 eS 9 10 | 11t 
m a b é h(d) M3 fs h(ds) h(d; f;*) Q Ror Ps 
(Series A: m= 1 (mod 4), ¢ = 1,d = 4m, d; = 4m;.) 
2 1 1 |-1| 1 | -6| 3 | 2 6 | -- 41 
3 2 1 }+1| 1 | -1/9 1 6 2 3 
6 5 2/+1/} 1] -2/9 1 6 2) -3 
7 8 3 |4+1/ 1 | -21] 3 4) 12 | -2 7 
10 o; 1) <2) 2) —mi 2 | 4] 12 | --- 157 
| | 
| | | | | 
11 0; 3 |+41| 1 —44| 3 4| 12 2 11 
14 6) 4/41] 1 —56 | 3 4} 12 | -2 7 
15 4; 1|+1/ 2] -5| 9 2] 12 3*| —5* 
19} 170| 39 | +1] 1 | —57] 3 4 | 24 2 19 
22; 197; 42 | +1] 1 | —66| 3 8 | 24 2) -11 
| | 
23 24 5 | +1) 1 —69 | 3 8| 2% | -2 23 
26 | 5 1 |-1/ 2 | —78} 3 4] 12 7 37 
30; IW} 2 |+1] 2 | -10| 9 2) 2% 5| —6 
31 | 1,520 | 273 | +1) 1 —93 | 3 4 12 -2 31 
34 35 6 | +1| 2 | —102) 3 4| 12 | -2 17 
35 | 6 1 | +1] 2 | —105| 3 8 | 24 5*| —7* 
38| 37 6 | +1] 1 | —114| 3 8 | 2% 2| —19 
39| 25 4/41] 2 | -13| 9 2) 2 3 | —13 
42} 13 2|+1|} 2 | -14|] 9 4 | 24 6| —7 
(Series B: m = 1 (mod 4), c = 2,d = m, ds = m3.) 
5 | 1 1 |-1]| 1 —15| 3 2 6 61 
13 | 3 1 |-1]| 1 —39 | 3 4| 12 29 
17 | 8 2|-1] 1 —51 | 3 2 6 13 
21 | 5 1 |+1/] 1 —7/ 9 1 12 3| -7 
29 | 5 1 |—1] 1 —87 | 1 (iii) 6 6 eee 13 
33| 46 8 |+1} 1 | -11| 9 l 6 | -3 11 
7| 8 2 |-1| 1 | -111] 3 8 | 24 . 53 
41; 64/; 10 | -1| 1 | —123] 3 2 6 | --- 5 
53 | 7 1 |-1]| 1 | —159| 3 10 | 30 | --- 17 
57| 302) 40 | +1/ 1 —19| 9 I 12 3 | —19 
61| 39 5 |—-1]| 1 | —183| 3 8 | 24 59 
65 | 16 2 |-1| 2 | —195!| 3 4| 12 29 
69 25 3 |+1]| 1 —23 | 1 (i) 3 3 | -3 23 
73 | 2,136 | 250 | -1| 1 | —219| 3 4| 12 37 
77 9 1 | +1] 1 | —231| 1@i) | 12/ 12 7 | —-11 
85 9 1 |} —-1] 2 | —255| 1q@i) | 12] 12 | --- 101 
89 | 1,000, 106  —1/ 1 | —267| 3 2 6 | --- 73 
93 29 3 | +1] 1 —31/ 1 3 3 3 | —31 
97 (11,208 1,188 | -1| 1 | —291) 3@) 4) 12 53 





* Here S = 2. (See Section 6). is 
+ When e = —1, Column 11 has the earliest prime p, for which H(p,) = 4. 
(See Section 6). 
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TaBie II 
Some Special Values of p for Which n | H(p) 

The table gives the minimum odd prime p(<1000) for which H(p) 
H(p) = 2n, if n is odd and e = N(e) = +1). If no such p occurs, the table lists 
pr the earliest prime (< 1000) for which H(p)/n (or H(p)/2n) gives the minimum 


n, (or 






































quotient r. 
| 
m n=8 16 | 32 9 27 5 | 23 | 7 | 11 
Series A 

2 137s | -353 | 269, | 79 | 643 | 199 
3* 313 193 181 | --- | 71 | J see | eee 
6* 409 97 | --- | 89 | 971, | 311 | | 743 | 109 
7° 71 751 | 127 | 179 | 271 | 131 | 197 | 617, 
10 241 449 ++» | 271 | -- | 19 | 419 | 131 

| Zee a 
11* 97 | 881 | 449 | 719 | 409 | 199 | 421 | 
14* 71 | 79 | +++ | 251 | ss | 29] + | 97 | 
15* 31 |: 163 | 487 | 61 | | 71 | -:: 
19* 73 --- | --- | 991 | 269 | 31 | --- | 13 | 397 
22* 353 401 | 641 | 883 | 593 | 271 | 701 | 127 | 131, 

| | | 
23* 4l 47 | 521 | | 59 | | 631 | 
26 641 881 oes! 139 | ++ | 3372 | 
30* 23 383 739 | 439 | 349 | 211 | 
31* 7 193 883; | 19 | 449 | 13 | --- 
34* 23 911; 163 | 59 | -++ | 83 | 433, 
35* 47 449 ms 71 | --+ | 89 | | 701 | -- 
38* 137. | 769 | 37 | 701 | 431 | | 127 | -- 
39* 673. | 79 | $27 | --- | 151 911 | 857 
42* 103 |: 673; 809 | 431 | 491 | 433 | -- 

Series B 
5 89 | se | 919 | 211 | 307 | 967 
13 233—C || :=--- | 827 59 | 211 | 109 
17 281 tee | 127 79 | vee | eee 
21* 199 337 | «- | 101 | 433 | 263 
29 233 =|: 673 971 | 619 | | 601. | 461¢ 
| | 

33* 71 =| 47 | «+++ | 43832 | 379 | 139 | | 239 | 331 
37 as | ee sss | TB | ses | ZL | oe are 
41 769: | 769 | --- | 307 | | 199 | | 491 | 593: 
53 929, | 929 | 449 | 433, | | 379 | 1132 | 659 
57* 487 | 127 | --- | 197 | 271 | 43 | -:- 
61 937 | 977 | --- | 271 | 487 | 59 | 463 | --- 
65 601 --» | 353 | 467 | 431 | 211 | eet or 
69* 71 | 230 | --- | 307 | --- | 79 | 97 | --- 
73 857 oe | | 107 | 379 | 3332 | 67 
77* 127 113 | vee | 101 71 | «+: 
85 ve vee | osee | TI 331 139 | 947 
89 809 | 641 | 929 | 631 | --- 59 | 503 | 967 
93* 463 | 79 | --- | 379 | 811s | 251 29 | 947 
97 1132 | 113 | 673 | 107 | --- | 151 — 


463 





(* Denotes values of m for which N(e) = 
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Tasie III 
Quadratic Forms Which are Perfect Cubes 
These are the forms (A, B,C) of discriminant B? — 4AC = d, f? which represent 
those primes p(+ 6m) for which 3 | H (p), where p is “conveniently” small. 
































m | ase | 4a | B| p |e] 4] B | p | ww 
Series A 

2 | —216 1 0; 7; 3); 2 or m1 fe 
3 | —324 1 0| 97) 12 | 2 2| 41 6 
6 | —) 1 0| 163; 6 | 2 0; 8 | 12 
7 | —756| 1 0; 193/ 12 | 7 0| 139) 6 
| | 2 2} 107| 6 | 14] 14| 17| 6 

10 — 1080 1 0; 271; 9 | 2 0| 137| 6 
5 0| 59 | 3 | 10 0| 37 12 

ll —1188| 1 o| 313! 24 | 0| 71 6 
| | 2 2| 149] 6 19 16| 19 6 

14 —1512 1 0| 379| 42 2 0; 191 | 24 
7 0; 61| 6 | 14 o| 41/ 6 

15 —1620| 1 0| 409) 12 | 5 | O| 1] 6 
2 2| 227| 12 | 10 10| 43 6 

19 | —2052 l 0| 577| 36 | 19 0| 103 6 
| Se 2| 257| 6 | 23 8| 2] 6 

22 —2376 1 0| 619; 6 | 2 0| 347 | 12 
il Q| 227) 12 | 22 | 0| 331 6 

| 7 +2 7) oe 14 | +12| 47 6 

23 | —2484/ 1 0;| 877 6 23 0} 131 6 
2 2| 311) 24 | 25 4|\ 349 6 

| 5 +4 5| 6] 10 +6 | 67 6 

26 | —2808 I >| wi 3 2 0| 353 6 
13 0 | 67 3 26 0 53 6 

30 | —3240 l 0| 811 | 30 2 0 503 12 
5 0| 167| 12 | 10 0| 241 60 

| 11 +4] 11 6 | 22 +4| 37 6 

31 | —3348 l 0} 853| 6 | 27 0| 139 6 
2 2; 419| 30 29 4 29 6 

34 | —3672 1 o| 919| 6 | 2 0| 461 | 6 
17 0| 71) 24 | 27 0! 61 6 

35 —3780 1 0| 1009| 12 5 0| 269 6 
7 0} 163} 6 | 27 0| 167 12 

2 2| 557| 6 | 31 8| 31 6 

10 10| 97| 6 14 14| 71 18 

38 —4104 1 0 | 1051 si? 0| 521 6 
19 0| 73| 2 | 27 0| 179 36 

23 +6 23 6 | 3 +22 31 6 

39 —4212 1 0| 1069 | 12 13 0| 337 . 
2 2| 5871 6 | 2% 26 | 47 6 

17 +2} 17] 6 | 31 +2/ 31 6 

42 —4536 1 0| 1303 | 6 2 0| 569 6 
7 0| 337| 8 | 14 0| 137 6 

13 
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TaBie II1]—Continued 
| j 
m | ds fz | A B p | H(p)| A B | p | H(p) F 
primy 
Series B 
5 | —135| 1 1 139 3 5 5 | 47 3 
13 —351 1 1 367 3 10 10 | 79 3 dita 
8 +1 11 3 A Breen “ 
17 —459 1 1 127 9 11 5 11 3 
21 — 567 1 1 571 6 7 7 | 109 12 
8 +3 23 6 vas | een ‘s 
29 —87 1 1 103 3 3 3| 41 6 
33 —891 1 1 223 6 1] 11 23 12 
37 —999 1 1 | 1063 3 16 5 | 619 3 
2 +1 131 3 4 +3/ 73) 18 
8 +5 89 6 SP Be F * see 
41 —1107| 1 1| 277| 12 | 17 7| 71 3 
53 — 1431 1 1 | 1447 3 20 13 | 239 3 
rf +5 7 3 18 +15 | 23 3 
| 10 +3 43 3 8 | +3 | 83 3 
57 — 1539 1 1 397 36 19 | 19} 1389 6 
5 +1 5 ie ee ee 
61 — 1647 1 1 | 1663 3 | 22 17 | 271 9 
18 +3 23 - kee +7} 53 j 
| 43 +11 13 6 e ee 
65 —1755 1 1 439 7 i @ 19 | 23 3 
5 5| 89| 18 | 13 13| 37 6 
69 —23 1 1 101; 6 | a Ga 
73 —1971 1 1 499 a | 2 23 79 3 
5 +3 5 6 tee 
77 —231 | 1 1 331 6 | 3 3 89 18 
8 5| 233) 6 | 7 7| 61 6 
85 —255| 1 1} 271| 15 | 8 1| s3| 21 
ii. 3 97| 12 | 5 5 | 131 3 
89 — 2403 | 1 1 601 | 12 24. 27 83 3 
93 mani 1 47| 6 ia rine es zi 
97 | —2619 | 1 | 1; 661; 12 | 27 27 | 3il 3 = 
3: | 27| Bi 3 noe va We ae 
bef 
and when m = 89, this holds for Dey 
Uni 
p = 53, 73, 109, 157, 233, 257, 269, 449, 461, 509, 601, 613, 641, 733, 757, 809, Tuc 
821, 929, 937, 977. | 
| TE 
Curiously enough, when m = 37 all p( <1000) for which H(p) = 0 (mod 4) ae 
satisfy H(p) = 4; from Table II, this value of m seems most ‘‘resistant to variety” _ 
in the values of H(p). out 
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TaBLe IV 
Irregular (Odd) Primes < 1000 


For values of m in Table I. Primes cf index 2 are marked with (*), unmarked 
primes are of index 1. (See Section 8.) 

















| | 
m p | H(p) m p | H(p) 
i 
Series A Series B 
2 13 2 13 241 2 
2 31 1 29 3* 1 
3 103 2 29 | ll l 
6 3 1 33 3 i 
6 7 2 33 29 2 
| | 
10 191 5 33 37 4 
10 643 1 37 7 I 
15 3 1 37 89 6 
15 181 2 37 257 6 
19 79 2 41 29* 2 
| 
oo Ge coy 7 oo ae 53 2 
22 73 2 | 53 5 2 
23 7 2 | 57 59 2 
23 733 2 | 69 5 2 
31 157 2 | 69 17* 2 
34 37 2 | 73 5* 6 
34 547 26 73 7 I 
35 23 2 73 4] 2 
38 5 2 ! 85 1 
39 5 2 | 89 5* 2 
| i 
39 7 2 89 7 1 
39 37 2 | 89 13 2 
42 | 3* Banorl 89 59 5 
42 | 5 2 | 93 13 2 
42 43 2 | 97 17 2 
| y 
42 | 71 2 





It is our hope that additional motivation might be suggested by these data 
before the next electronic tour de force is attempted. 
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A Very High-Speed Digital Number Sieve 


By D. G. Cantor, G. Estrin, A. S. Fraenkel, and R. Turn 
1. Introduction. The general sieve problem may be stated as follows {3}. Let 


m, M:,,---, m, be s positive integers, relatively prime in pairs. Consider the 
congruences 
(1) x = a;;(mod m,), 62 1,2;++- ,e37 = 1,2, -** 5 < m;. 


For fixed i, the a;; are distinct non-negative integers less than m;. The problem 
is to find all integers N between given limits, say 


(2) ASN <B, 


such that N is a solution to s of the congruences. (It is, of course, clear that no N 
can be a solution to more than s of the congruences (1).) 

Examples: On the one extreme there is the Sieve of Eratosthenes for finding all 
primes p in the range A = B’* < p < B, where t; = m; — 1 for all i. (Here m; 
are all the primes < B’’*.) On the other extreme there is the Chinese remainder 
type of problem, where ¢; = 1 for all i, and there is only one solution among [| ]j.1 m, 
numbers. 

In between these two extremes, there is the important quadratic sieve, 
where roughly ¢; = m;/2 for all 7. It is used in problems involving quadratic residues, 
Diophantine equations of second degree and other quadratic type problems. 

About thirty years ago, Lehmer [1], [2] constructed a novel special-purpose 
device for sifting. It used the first 30 primes as moduli. Its processing rate was 


3 X 10° numbers/min. 


General-purpose computers are not very well suited to sifting, and the earlier 
models could not compete with Lehmer’s machine. However, the speed of the more 
recent machines makes up for their lack of orientation towards the sieve problem 
insofar as surpassing the performance of Lehmer’s machine is concerned. Thus, 
the rate for a quadratic sieve using the first 30 primes on the IBM 7090 is approxi- 
mately 


7 . 
10° numbers/min. 
The present paper describes a special-purpose device, where rates in excess of 
10 . 
10° numbers/min. 


can be achieved for quadratic sieves. The device consists of basic digital building 
blocks from which a suitable sieve is assembled for each problem. Thus, by an 
appropriate rearrangement of the building blocks, problems with different moduli 
can be run. It is also shown that if the device contains a certain minimum amount 
of hardware and is attached to a general-purpose computer, then problems can be 
run where, roughly speaking, the number of moduli is not limited any more by the 


Received June 23, 1961. The preparation of this paper was sponsored, in part, by the 
Office of Naval Research and the Atomic Energy Commission. 
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amount of hardware of the special-purpose device, but only by the size of the 
memory of the general-purpose computer, and the rate is still of the order of 10° 
numbers/min. We use a so-called “Fixed Plus Variable Structure Computer” 
organization for realizing the combination between the special- and general-purpose 
computers [6]. This also enables one to use the digital building blocks of the sieve 
for building other special-purpose devices which one might want to associate with 
the general-purpose computer. 


2. Binary Set-Up of the Sieve. For solving the system (1) on a digital computer, 
we consider a matrix M of sizes X (B — A) withentriesc,, (i = 1,2,---,s;k =A, 
A +41,---,B — 1), defined by 


lifk= Qi; 
Ce = 
0 otherwise 


(mod m;) 
(j = 1,2, --- ,&). 


Then every column, all of whose entries are 1, corresponds to a solution, and con- 
versely. 


Example: Find the primes p such that 
(3) 6 Sp < 36. 


The relevant congruences are x = 1 (mod 2), x = 1, 2 (mod 3), x = 1, 2, 3,4 
(mod 5). The matrix M is given by 





XK i. ep ages: Fe — 





for Ny " 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 

. 

| 

| NI ah. = ‘3 
2 P22 SF ECP OBA ARARRHCOA@1O1@ tet O11 Oe 1 
3 aw 2s es eee eee eee rer 21 eakernenka 
5 Pees Cie © ee See eee ak OcT Aaa Ok 841t 1 @ 





The columns all of whose entries are 1 correspond to the primes in the range (3). 

The rows of M are periodic with period m;. Thus, the first ordered m; bits of 
the ith row determine the rest of this row completely, and we call them the periodic 
pattern e; of m; . 


3. Method of Solution on the Special-Purpose Computer. We now give an 
informal introduction to the principle of operation of the special-purpose device. 
It will be observed that the method is based on ideas used in earlier work [1], [2], 
[3] in this field. 

A first approach to the mechanization of a special-purpose sieve would be to 
build a matrix precisely in the form displayed for the example above with observa- 
tion posts in every column, detecting coincidences of non-zero bits. Problems which 
can be solved by such a procedure are limited to those which can fit into the maxi- 
mum size matrix which can be assembled, i.e., computationally trivial problems. 

We order the moduli so that 


(4) m <M <-*+ mM. 
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Since solutions exist only corresponding to columns with non-zero bits, we may 
eliminate the m,-row and many of the components required to detect coincidences, 
by establishing coincidence gates only in those columns where the m,-row has 
non-zero entries. 

Large problems may be handled by constructing only m, columns of the matrix 
and testing for solutions in these columns in parallel. Entering the next batch of m, 
numbers turns out to be equivalent to performing prearranged circular shifts in 
the s — 1 rows. This procedure would require a matrix of size (s — 1) XK m, with 
coincidence gates established in the columns as prescribed above. It is possible to 
use such a matrix to define potential solutions even when it is only feasible to 
mechanize | < s — 1 rows of the matrix, and then have a general-purpose computer 
complete the test for solution. 

The range of problems which may be handled is increased when it is recognized 
that the periodic pattern e; associated with the ith row completely determines the 
rest of the row. In the following we give an algorithm which defines a procedure 
requiring only m; elements in the ith row, giving up only the regularity of the 
coincident gate connections. The special-purpose computer consists of basic digital 
building blocks or modules which are assembled into a matrix consisting of s — 1 
shifting registers, the ith of length m; , and initially containing the periodic pattern 
e: (t = 1, 2,--- , s — 1). Observation posts are placed at certain positions in the 
matrix which sift out the solutions to (1) among the first m, numbers.* Next, a 
circular shift is performed in each register, which is equivalent to bringing in the 
next m, numbers to be sifted. This is followed by the observation posts sifting 
out the solutions among this new batch of numbers. This process of sifting followed 
by shifting is continued until all the numbers are processed. 


4. The Algorithm. We divide the numbers N in the range (2) into sets S,, defined 


byt 
‘. . B-A-l1 
(5) S,={N:N<B,N=A+nm,+k;0sk<m,}, n=0,1,---, ——e 
Thus, each set (except possibly the last) contains m, numbers. 
Let 
(6) m= Qm+r, 0O<r<m (*#=1,2,---,8s—1). 
With each set S, we associate a matrix M, of size (s — 1) X m, with 
entries c;;(n) (¢ = 1,2,---,s — 1;j = 0,1,--- , m, — 1), defined recursively by 
(lifO<j<m and A+j =ai;,---,a:,4 (mod m,) 


(7) (0) = 4 
\0 otherwise. 


(c:jar(n —1) if OSG +r: < m 
(8) Cis(N) = 4Cijtr—m(n — 1) if OSZ<m and j+r42= mM, 
(0 if mSj<m,. 





* The positioning of the observation posts is determined by m, and its residues in such a 
way that a register of length m, is not required. 
Tt [x] stands for the largest integer < z. 
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Equation (8) can be written in the form 
ca(n — 1) where m >d=j+r;(modm,), if 057 < m; 
aie if m Sj <m. 
Hence (7) and (8) are equivalent to 
\ f# Ofj7<m; and A+ j =a4;,., — nr; (mod m) 
(9) c:j(n) = 
0 otherwise (9, = 1,2, --> , &). 
If N « S, is a solution to the system (1), then by (5), 
(10) A+k=a,,.,(modm,) (0 £k < m,;v, = 1,2, ---, t,) 


for all n. 
By (5) and (6) we have also 


(11) A+k=a,;,, —mr;(modm;) (vv; = 1,2,---,t#;¢ = 1,2,---,s— 1). 
Hence, if we let 
(12) k= wm: + ui, 0S ui <m; (4 = 1,2,---,s— 1), 
then by (11), 
A + ui = Gi», — nr; (mod m,) (v; = 1,2, ,t;i = 1,2,---,8— 1), 
so that 
(13) Cru (nm) = 1 
fori = 1,---,s — 1 by (9). 
Also the converse holds. That is to say, if (13) holds fori = 1,---,s — 1 
(where u; is given by (12) and k by (10)), then 
(14) N=A+nm+t+k 
is a solution to the system (1). This is the basis of the algorithm. We list the ¢, 
solutions k, , --- , k:, of (10), and for each of them its corresponding s — 1 values 
u;. Then the numbers N = A + nm, + k,, (v, = 1,--- , t,) for which (13) holds 
fori = 1, --- ,s — 1, are solutions to (1), and these are all the solutions. 
Example: Find all solutions in the range 
—-15 SN < 25 
to the system of congruences 
z=1 (mod 2) 
z=1,2 (mod 3) 
x = 2,3,4 (mod 5) 
xz =0,1,2,5 (mod 7) 
z=0,1,8,9 (mod 11). 


Reference is made to Table I. The matrix M, is constructed according to (7) 
(omitting all strings of zeros). Equation (8) implies that the 7th row of matrix 
M,, is obtained from the ith row of M,_, by means of circularly left shifting the 
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TABLE I 

The Matrices for the Problem 
ag | 18 <6 <8 8 ar =90" 9 387] be os 
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first m; bits by r; positions (or by circularly right shifting them by m; — r; posi- 
tions). In the present case, r; = 1, re = 2,73 = 1, 7% = 4, so that in passing from 
one matrix to the next, the first row is shifted left circularly by 1 position, the 
second by 2, the third by 1 and the fourth by 4 positions. This is the way M, , 
M, and M; are obtained. 

The values k, , ke, kz, ks of Table I] are computed by (10), and the corre- 
sponding values of u; by (12). Table II defines a pattern of observation stations 
which sift out the solutions in each matrix; the entry u; represents the column co- 
ordinate corresponding to the row coordinate i, at which an observation station 
exists. In My, the observation pattern u; = 0, 2, 2, 2 indicates a solution, since all 
matrix positions corresponding to that pattern are filled by 1’s. They appear in 
bold type in Table I. The corresponding value of k is k, = 2. The solution is, there- 
fore, N = —15 + 0 X 11 + 2 = —13 by (14). The only two other solutions 
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Fig. 1—Conventional Shift Register. 











are found in M;, for kj = 1 and kz = 5, also indicated by bold type in Table I. 
They areN = —15+3 X 11+1= 19andN = —15+3 X 114+ 5 = 23. 


5. The Special-Purpose Computer. The special-purpose device should be 
flexible enough to allow usage of different moduli. Therefore the basic building 
blocks of the device are modules consisting of memory elements and gates which 
will be assembled into the appropriate-size shifting registers and sets of and-gates 
for each problem. The and-gates are the realization of appropriate combinations 
of the observatior posts, and test for coincidence of 1’s. 

The example of the previous section suggests the construction of the device. 
Its central part consists of shift-registers R,,--- , R,., the ith of length m; , 
which will store the matrices M,, (without the trailing strings of 0’s). Register R; 
will shift circularly left by r; positions. It is important to note, however, that this 
can be effected in one shift time, rather than in r; shift times, and further, that the 
wiring can be so arranged that any transmitting memory element is adjacent to 
its receiver. In order to do this, we rename the memory elements in R; so that 
element number 0 is at the left, followed by element number r; , followed by num- 
ber 2r; (mod m,), by 3r; (mod m;), etc. Then each transmitting element is adjacent 
to its receiver, and every element will appear exactly once. Fort (m;, m,) = 1, 
so that also (m;, 7:1) = 1 (¢ = 1,---,8s — 1). Hence 7; generates the additive 
cyclic group of non-negative integers (mod m;) and every element appears exactly 
once in the register. 

Ezample: Suppose that for a certain sieve problem m, = 23, and m; = 10 for 
some i < s. Then R; has to shift circularly left by r; = 3 positions. On a conven- 
tional shift-register, three shifts of the type indicated in Figure 1 would have to be 
performed. The same result can be obtained in one shift time by the specially 
wired-up register of Figure 2. However, the long wiring has an undesirable effect 
on the speed of the system. Renaming the memory elements as in Figure 3, the 
three shifts can be done in one shift time with conventional wiring. 


t (a, b) stands for the greatest common divisor of a and b. 
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Fig. 2.—Specially Wired-Up Register for Performing a Shift of Three Positions in 
One Shift Time. 
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Fig. 3.—Final Form of the Register. 


















































The second part of the special-purpose device consists of t, sets of and-gates, 
one set for each value of k which is a solution to (10), “anding”’ together positions 
c;; in the registers, as defined in the previous section. Also a counting register R 
with capacity > (B — A)/m, and a shift control are required. (The shift control 
does not normally have to be reconstructed for every problem.) The register R 
contains the value n of the matrix M,, currently being sifted. 

The sifting process consists of the following steps. 

1. Clear R. 

2. Load R,, --- , R, with Mz , whose entries are defined by (7). 

3. Record n and record the values k for which coincidence is obtained, i.e., 
the values k associated with the sets of and-gates which are excited, if any. The 
corresponding solutions are given by (14). 

4. Advance R by unity and perform a circular shift in each shift-register 
according to the above scheme. This effectively reloads the registers with the next 
batch of m, numbers to be sifted. 

5. Terminate process ifn > (B — A)/m, . Otherwise go back to 3. 

It is thus seen that in this process m, numbers are processed per shift time. 


6. The Sieve in a Fixed Plus Variable Structure Computer. It was remarked 
by Lehmer that special-purpose equipment attached to the arithmetic unit of a 
fast computer can speed up computation of permutation problems [4], and of other 
problems [5]. More generally, we consider a so-called “‘Fixed Plus Variable Struc- 
ture Computer” (to be designated by (F + V) computer), which consists of a 
conventional digital computer (the fixed part to be denoted by F), and a set of 
modules (the variable part to be denoted by V). Many problems contain a part 
which can be solved on a special-purpose computer in a much more efficient way 
than on a general-purpose computer. For such a problem, the modules are assembled 
into a suitable special-purpose device which handles this part. The rest of the prob- 
lem is handled by F. A supervisory control coordinates the operation of the two 
computers. However, the special-purpose configuration is not retained permanently, 
but may be reorganized into other configurations for other problems. For a more 
detailed description of the concept of the (F + V) computer, the reader is referred 
to the literature [6]. 

The special-purpose device described above has a limited amount of hardware. 








148 D. G. CANTOR ET AL 




























































































V F 
J 
————— __ a > 
Rx, o } ny epee 
= cs 
SHIFT SUPERVISORY F—® CONTROL Re 
CONTROL | -+ e+ CONTROL - 
_ MEMORY 
@s 
_ —$——————— ) 
Re- 
k- REGISTER 
ARITHMETIC 
OR 






































Fic. 4.—(F + V) Computer Organization for the Sieve Problem. 


For certain problems it may be desirable to use more moduli than can be mech- 
anized with the available hardware. In order to handle such problems, we imbed 
the special-purpose device in an (F + V) computer. This allows, as will be shown 
subsequently, handling problems in which the number of moduli is limited only 
by the number of periodic patterns e; of m; that can be stored in the memory of 
F, provided that the hardware of V is sufficient to mechanize the first | of the s 
moduli. The parameter / depends on the relative number of 1’s and 0’s in the 
periodic patterns of the sieve, and on the relative speeds of F and V. The use of 
F + V also allows using the modules for building special-purpose devices other 
than the sieve, and attaching them to F. 

Figure 4 shows the organization of the (F + V) computer for the sieve problem 
by means of a block diagram. The V-part, which acts as the special-purpose device, 
mechanizes the first / moduli. The k-register records the values of k (the solutions 
of (10)) corresponding to coincidence. The periodic patterns e;4,;,---, e of 
Mi41,°** , m,are stored in the memory of F. Numbers will be sifted in V, and when 
coincidence occurs, the contents of R and of the k-register are transferred to F, 
where so-called solution candidates N of the form (14) are formed, one for each 
value of k. Then divisions of the form 


(15) N = om; + pi, 0S pi < m; (¢=1+1,---,8) 


are performed in F. The residue p; determines uniquely the position of the bit of 
the periodic pattern e; of m; corresponding to N. The number N is a solution if 
and only if these bits are 1 fori = 1 + 1,--- , s. Thus, the (F + V) computer 
is so organized that V will do the high-speed sifting, and F will do divisions. | will 
be chosen so that the average time per coincidence in V is at least as large as the 
average division time per coincidence in F. Then V would normally do its divisions, 
until such time when coincidence is obtained in V. At such time, V is interrupted 
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and the transfers from V to the memory of F occur, the latter acting as a buffer, 
capable of storing “bursts” of solution candidates which V might produce occa- 
sionally. After the transfers, both parts are again decoupled and assume their 
respective tasks. The program for V is outlined in the following six steps. 

1. Clear R. 

2. Load R,,--- , Ris with My. 

3. Check for coincidence. If none is obtained, go to 5. Otherwise continue with 4. 

4. Interrupt F and V. Transfer the contents of R and of the k-registers to the 
memory of F. 

5. Advance R by unity and perform a circular left shift of r; places in R; (i = 
1,---,l—1). 

6. Terminate process ifn > (B — A)/m,. Otherwise go back to 3. 

The program for F is simply to produce solution candidates N of the form (14), 
and to perform divisions of the type (15) for each of them until the first 0-bit is 
encountered. If none is encountered, N is recorded as a solution. 


7. Speed and Hardware. The speed and hardware requirements will now be 
discussed in terms of an example, for which we choose a quadratic sieve problem 
where the moduli are the first s primes. The first column of Table III contains 
values of the independent variable |, the number of moduli mechanized in the 
special-purpose device. The table displays the speed and hardware requirements 
for such a sieve as a function of |. The second column contains the /th prime. Let 
t be the total time required for performing the coincidence test and the subsequent 
circular shifts in the registers. Using the register organization described in Section 
5, the circular shifting amounts to a left shift of one position in each register. We 
assume transistorized circuitry, for which 


t = 0.2 uw sec 


is chosen (speeds approximating those of the IBM 7090). Thus, m; numbers are 
processed in this time if no coincidence occurs. (If the registers are of the double- 
rank type, both ranks will be equipped with sets of and-gates, and ¢ = 0.2 u sec 
is the time for a coincidence check and for transferring one rank into the other. 
Thus also in this case m; numbers are processed in 0.2 » sec.) 

We consider first the case s = 1, that is, we use only a special-purpose computer 
without a conventional general-purpose computer. Assuming the solutions to be 
sparse, so that we may neglect the time of recording them, the rate of the sieve is 


v = 6 X 10’ X m, x mi 


; = 3 X 10° X m, numbers/min. 


These values are displayed in the third column. 

Since the sieve is quadratic, the probability of any randomly selected bit in the 
periodic pattern e; to be 1 is about 0.5. Hence, on the average, one coincidence is 
obtained per 2’ numbers sifted, or every 


These values appear in the fourth column. 
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If s > l, divisions have to be performed in F. The probability that exactly i 
divisions suffice to decide whether any solution candidate N has to be rejected or 
accepted is (14)* (1 < i < s — 1). Hence the expectation of the number of divi- 
sions for each N is given by 


s—l 

(16) y= 2d, 1/2' = 2—(s—1+2)/2°". 

Thus the average number of divisions for each solution candidate approaches 2 
asymptotically from below. Assuming the IBM 7090 as the fixed machine F, this 
division subroutine takes about 200 u sec for two divisions. Also, preliminary studies 
of the mode of transfer from V to the 7090 indicate that the transfer of the con- 
tents of R and of the k-register requires no more than 7 yu sec. (See appendix.) That 
is, this is the maximum time during which V is idle. F is interrupted only insofar as 
it requires memory access during this time. Actually, V could already resume its 
operation after the transfer of n from the R-register. It would have to wait additional 
time only if a new solution candidate is formed before the current contents of the 
k-register has been stored away, which is a rare event. However, in our computation 
of the overall speed of the sieve we assumed that V is interrupted for 7 u sec. during 
each transfer. 

Thus, the average overall rate of the sieve is given by 


v ; 
w= iti numbers/min. 

for r 2 200 uw sec. If r < 200 uw sec, V will have to wait for F, and the average 
overall rate for this case is 


T v 


~ 200 1+ 7/r 


Thus, the operation of the sieve becomes rapidly more and more inefficient as 7 
decreases below the critical value of 200 u sec. The values of w appear in the fifth 
column of Table III. The last column displays the required number h = >< {=i m, 
of memory elements for the special-purpose device. (This number has to be doubled 
if the registers are of the double-rank type.) 

The lower bound for / in Table III was chosen to be 9 because for 1 = 8 the 
rate would already be less than can be achieved with conventional present-day 
computers. The upper bound was chosen by setting arbitrarily a hardware con- 
straint of 1500 memory elements. 

The lowest value of 7 for which r = 200 uw sec is r = 247.3 uw sec. Thus V should 
contain at least 15 registers consisting of 328 memory elements. Figure 5 displays 
the overall rate as a function of required memory elements. Two simple conclusions 
can be drawn from the monotonicity of w as displayed in Figure 5. First, / should 
be chosen as large as possible. That is to say, as much hardware as available should 
be thrown in to build the sieve; even so the cooperation between F and V becomes 
less efficient as / increases beyond the critical value of 16, in the sense that F be- 
comes more idle. Secondly, the use of slower memory elements is indicated if a 
larger number of them is available, hence the possibility of using magnetic core 
registers. 


w numbers/min. 
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Tasie III 





Speed and Hardware as a Function of | For a Quadratic Sieve 




















I—No. of Moduli . | h—Number of 

Implemented in — 4 meg | r Average Time | Lg ie | | Memory Ele- 

" tae) Modulus | Sieveifs =} | Pe Coincidence | Sieveifs>1 ante Thevice 

| 

9 | 23 6.9 xX 10° 4.5 uw sec 5.70 X 107 | 77 
10 | 29 |8.7 xX 10° 7.lusec | 1.56 X 10 | 100 
11 31 9.3 xX 10° 13.2 uw sec 4.03 X 108 129 
12 a 1.11 X 10° 22.1 uw sec 9.28 x 108 160 
13 = 1.23 X 10'° 39.9 uw sec 2.09 X 10° 197 
14 | 43 |1.29X10%| 76.2usec | 4.49 x 10° 238 
15 = 1.41 XK 10'°| 139.4 uw sec 9.34 X 10° | 281 
16 53 1.59 K 10° | 247.3 uw sec 1.54 X& 10° 328 
17 | 59 1.77 X 10" | 444.3 uw sec 1.73 X 10° 381 
18 | 61 |1.83X10"| 859 ywsee | 1.81 X 10” 440 
19 | 67 2.01 X 10%° 1.6 msec | 2.01 X 10% 501 
20 71 2.13 X 10" 2.9msec | 2.13 X 10 568 
21 | @ 2.19 X 10% 5.7msec | 2.19 X 10% 639 
22 i; 79 2.37 X 10% 10.6 msec | 2.37 X 10' 712 
23 | 83 2.49 X 10% 20.2 msec | 2.49 X 10% 791 
24 | 89 2.67 X 10'° 37.7 msec | 2.67 X 10% 874 
25 97 |2.91 X 10°| 69.2msec | 2.91 X 10” 963 
26 101 3.03 X 10°; 132.9 m sec 3.03 XK 10° | 1060 
27 103 3.09 X 10% | 260.6 m sec 3.09 X 10% | 1161 
28 | 107 3.21 & 10° | 501.7 msec 3.21 XK 10'° | 1264 
29 | 109 3.27 XK 10° | 985.1 msec | 3.27 X 10! 1371 
30 113 3.39 XK 10'°| 1900 msec | 3.39 X 10° 1480 








By (16), the average number of divisions per solution candidate in F is less 
than 2, whatever the number of periodic patterns that are stored in F. Therefore, 
the rate w of a quadratic sieve is practically independent of s, and the number of 
the moduli is limited only by the number of periods that can be stored in F. A 
similar remark applies also for sieves that are “less than quadratic,” i.e., where 
the number of 1’s in e; is < m;,/2. For these types of sieves there are even less 
divisions to be performed, and a higher overall rate is obtained. For sieves that are 
“more than quadratic,’ and in particular for those which approach the type of 
sieve of Eratosthenes, more than two divisions are required on the average, and a 
higher critical value of / is obtained. 

Variations of the above described method which result in even higher speeds 
(and therefore involve higher critical values of 1) are clearly possible. For example, 
two moduli may be combined in V, say m, and m,_, , by initially solving the two 
congruences involving m; and m;_; manually or on F. Then mm;_, numbers can be 
processed per shift time, for which ¢;t;_, sets of and-gates are required. Registers m; 
and m;_, do not have to be built of course. As another example, we might mechanize 
the moduli m, , m2, --- , mi1, m, in V, rather than m,, m2, --- , mis, M;, SO 
that m, numbers instead of only m; are processed per shift time. 
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Fic. 5.—Rate of Sieve as a Function of the Number of Memory Elements. 


Similarly, two moduli may be combined in F. For example, combining m)4, 
with mi42 and m4; with mii, (which increases storage requirements in F), and 
using as divisors the moduli mji.mi42, Mi4sM144, Miss, ++, m,, the average 
number of divisions that have to be performed approaches 11/8 asymptotically 
from below. Thus, the effect of combining moduli in F is to lower the critical value 
of r. Such a procedure would therefore be used when the available hardware in V 
is smaller than required for keeping up with the speed of F implied by an average 
of two divisions per solution candidate. 

The high speeds which can be achieved by our method suggest its applicability 
for conversion of numbers from the modular number representation [7] to the 
conventional polyadic representation. Since this problem is of the Chinese remainder 
type, it seems possible to include in the sieve special solution hunting properties. 


8. Conclusion. A method has been presented to sift numbers satisfying a set of 
linear congruences from among a large set of numbers. The important properties 
of the resulting special-purpose device are that a relatively large set of numbers is 
processed essentially within the time required for performing a shift of one position 
in an ordinary shift-register, and that no memory references are necessary. This 
leads to an overall speed gain of about three orders of magnitude over modern 
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present-day computers such as the IBM 7090. By combining the device with a 
general-purpose computer, the size of problems that can be run is greatly increased 
with almost no decrease in speed. 


Appendix 
The Division Subroutine. For the purposes of this subroutine, written for the 
IBM 7090, we restrict the size of N in (2) to a number representable by 72 binary 
bits. 
The first 36 of these are called HW, and the last 36 are called LW. The core of 
the subroutine consists of the following sequence, where it is assumed that the 
accumulator is cleared at the beginning. Every bit of the periodic patterns e; is 


stored in a separate word, denoted by WM, and the corresponding period is stored 
in M. 


LDQ (Load the MQ) HW 
DVP (Divide) M 
LDQ (Load the MQ) LW 
DVP (Divide) M 


PAC (Place complement of address 0,4 

in index register) 
CLA _ (Clear add) WM,4 
TMI (Transfer on minus) OUT 


This sequence requires 36 cycles. For a quadratic sieve, the sequence has to be per- 
formed twice on the average for each solution candidate. Another 20 cycles are 
required for performing the multiplication and addition implied by (14) and 
bookkeeping. One cycle takes 2.18 » sec. Thus the subroutine requires about 200 
bu sec. 


Transfers from V to F. A preliminary study of the(F + V) organization based 
on the IBM 7090 as F indicates that transfers from V to F can be effected in the 
manner of a data channel. Such a channel has a “Channel Address Counter” 
(CAC), from which addresses are transferred to the ““Memory Address Register.” 

Suppose that the memory region bounded by addresses K and K + M is al- 
located for storing the value n contained in R, and L to L + M for storing the 
contents of the k-register. We assume, for simplicity, that registers R and k do 
not exceed 36 bits. In its normal form, the CAC contains an address of the form 
K +7i(0 si S M). Three flip-flops FF1, FF2, FF3 are contained in SC. FF1 re- 
cords whether F or V was the last user of the buffer region of the memory. FF2 
and FF3 define “full” and “empty” conditions of the buffer. 

We adopt the following operating rules. 

1. When V wants to store into the memory, the CAC is advanced by 1 if FF1 = 1, 
and remains unchanged if FF1 = 0. Then n is stored at the address currently held 
in CAC, say K + i. Next CAC is changed to L + i, and the contents of the k-reg- 
ister are stored. Then CAC is set back to K + 7. After execution of these stores, 
FF1 is set to 1. 

2. When F wants to fetch a pair of new values from the memory, the CAC is 
decreased by 1 if FF1 = 0, and is left unchanged if FF1 = 1. The address (con- 
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tents of CAC) is forced into the F Memory Address Register as a consequence of 
recognition of a special instruction in F by the Supervisory Control. Both the 
value n and the corresponding contents of the k-register are then fetched by the 
previously described K — L interchange, and CAC is set back to K + 7. At the end 
of the fetching operations, FF 1 is set to 0. 

Thus F always handles first the latest information brought in from V. If at any 
time CAC holds the address K + M, and if FF1 = 1, then FF2 is set, which pre- 
vents V from storing into the memory. (Of course for sufficiently large /, such an 
occurrence is very rare.) FF2 is reset by the resetting signal of FF1. If at any time 
CAC holds the address K, and if FF1 = 0, then FF3 is set, which prevents F from 
fetching. FF3 is reset by the setting signal of FF1. 

Preliminary studies of this mode of transfer indicate that transfer of the first 
word takes at most two cycles, and the second takes one cycle. Thus the transfer 
of n and the contents of the k-register from V to F requires approximately 7 u sec. 
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On Finite Difference Methods of Solution of the 
Transport Equation 


By R. P. Pearce and A. R. Mitchell 


1. Introduction. In recent years several difference schemes have been proposed 
for solving the transport equation 


ui, au _ 
(1) at V(z, t, u) = = F(z, t,u) 


in one form or another, where V is the velocity of propagation of a profile given 
initially along the z-axis. Most of these schemes can be found in Richtmyer [1] 
and, generally speaking, they are chosen primarily from the point of view of sta- 
bility. 

An equation of the type (1) has a single family of characteristics in the (rz, #) 
plane and in any step-by-step method of solution it is essential from the point of 
view of accuracy that the characteristics be followed as closely as possible. It is 
proposed to examine existing difference schemes from this standpoint and to derive 


new formulas of greater accuracy. For the purposes of this paper, it is sufficient to 
consider the simplified version of (1) 


ou ou 
(2) RAS dS 


where V is constant, from which it follows that the given profile at t = 0 is propa- 
gated without change of shape in the direction of the z-axis with velocity V. If a 
difference scheme fails to give an accurate solution of (2), it is pointless to con- 
sider it as a means of solving more complicated forms of (1), in particular, forms 
which incorporate variable velocity of propagation and source or sink terms. On 
the other hand, it is realized that schemes which successfully solve (2) may not 
give comparable accuracy when used to solve (1). In the case of (1), the char- 
acteristics are curved and can only be determined by integration of the equation 


= = V(z, t, u). In addition, the equation ~ = F(z, t, u) has to be solved. These 


computations, however, involve only numerical integration, a process which can 
be made as accurate as required in most problems. 


2. Stable Finite Difference Schemes Now in Use. Existing stable difference 
schemes will now be discussed with reference to equation (2). The characteristics 
of the latter are straight lines inclined to the ¢-axis at an angle 


(3) 6 = tan V. 
— V At 
In these schemes, the parameter p is introduced where p = as” and Az and At 
are the respective mesh lengths in the x and ¢ directions. 
Received January 10, 1961. 
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Fig. 1 


Difference System I (Friedrichs [2]). This is given by 
(4) Ur s+1 = 3(1 niet P)Ur+1,s + 3(1 + P)Ur-A,s 


where x = rAz and ¢t = sAt. This system can be obtained by replacing ~ and - 
Cc 


at the node (r, s) by sxe (tnt — Ur1,3) and Coe — Ur,s) respectively, then 


substituting }(U,+41,. + U1.) for u,,,. Another and more satisfactory way of de- 
riving (4) is now proposed. In Figure 1, the characteristic through P cuts AC in 
P, where BP; = pAz, and it follows that up = up, . Since P; is not a mesh point, 
the value of « at P; may be obtained by linear interpolation between A and C, 
and so (4) is obtained. In addition, since the coefficients on the right-hand side of 
(4) have sum unity, the solution computed by (4) is bounded if both coefficients 
are positive which leads immediately to the condition | p| < 1 for stability (Richt- 
myer [1], p. 43). 
Difference System II (Carlson [3]). This system is given by 


(5a) Urea = (1 — p)tr.s + Pirie (0sp<1l) 


be 1 p 
fae? to = = + 8, (p>1 
5b) +1 i+? , Il+p 1s+1 p ) 


and two similar formulas if p < 0. 
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In Figure 2, PP; is the characteristic through P. In this scheme, three points 
only are used, the choice of points depending on the position of P,. If 0 < p < 1, 
P, lies between B and C and the formula used is (5a), whereas if p > 1, P, lies 
outside BC and the formula used is (5b). It is presumed that these formulas were 


. eae . OU 1 ou 1 
ll b f — pers r,s pes r.s _— em i r— :) 
obtained originally by replacing 7 by XS +1 — Ur.) and = by An : C4. 


or (tron — Uri.) forO S p S 1 and p > 1 respectively. 


When P; lies between B and C (Figure 2a) it follows that BP,;:P,C = p:1 — p. 
Thus, on using linear interpolation between B and C together with the result 
Up = Up, , formula (5a) is obtained. When P, lies beyond C (Figure 2b) it can be 
shown that BP.:P,2R = p:1, and so using linear interpolation between B and R 
together with up = up, , formula (5b) is obtained. The solution computed by (5) 
is bounded for all p, as the right-hand sides of both (5a) and (5b) sum to unity 
and have positive coefficients. 

As Carlson’s scheme has been used extensively to solve problems involving the 
transport equation [1], [4], it is worth studying in some detail with a view to 
determining its probable accuracy. If 0 < p < 1, formula (5a) is as accurate as 
the linear interpolation of u between B and C. As these are neighboring mesh points 
on ¢t = sAt, the line of most recently computed values of u, it is to be expected that 
(5a) will give reasonably accurate values of u. Certainly (5a) will be superior to 
scheme (4) proposed by Friedrichs since the latter uses linear interpolation of u 
between A and C, mesh points two distance intervals apart. If p > 1, however, a 
much less satisfactory state of affairs exists. In Figure 2b, RR, is the characteristic 
through R, and theoretically ue = ue, . Similarly, wp = up, and (5b) becomes 
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] 
TP TF>5 wet pou Ri- 
Since BP,:P,R, = p:1, this formula is equivalent to linear interpolation between 
B and R, which for large values of p, where B and R, are many distance intervals 
apart, may be very inaccurate. In fact, the foregoing seems to suggest that implicit 
schemes in general are poor, particularly for large values of p. 
Difference System III (Central Difference Formula). This is given by 


(6) Ur s+ = Uris—1 — DUr+i.s + PUur-1,s 5 


and has been used with success by Malkus and Witt [5] to solve some problems in 
meteorology involving —" : temperature and oe 7 two dimensions. 


It is obtained by replacing = — by — Dal tas, e+1 — Ur,s—1) and S™ = by 55 i (te. a awe 


Alternatively, from Drea § 3, if GG, and PP, are the elena through G and 
P respectively, so that ue = ue, and up = up, , (6) may be written as 


Up, = Ue, — PU, + PUc. 


This result is equivalent to using a parabolic interpolation formula incorporating 
values of u at A, G,, and C, thus (6) is expected to be an accurate formula, par- 
ticularly for small values of | p |. In fact, (6) is stable for —1 < p < 1, and can 
only be used if G, and P; lie between A and C. 

It is interesting to compare the foregoing predictions of accuracy with numerical 
calculations carried out using difference systems I, II, and III in turn to solve (2). 
Two initial profiles of u are considered, the “roof top” and the “sine” and these are 


illustrated in Figure 4. All calculations are carried out until a time “os is reached. 
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Theoretical values are used at the second time step in order to start the calculation 
using System III. The results, accurate to 0.001, are shown in Tables 1(a) and 1(b) 
for the “roof top’”’ and “sine” profiles respectively. The last row of these tables 
gives the sum of the moduli of the errors }» | e |. The outstanding features of these 
results are the poor accuracy of Carlson’s scheme for | p| > 1, and the compara- 
tively high accuracy of the central difference formula. 


3. Two-Level Interpolation Schemes. As a consequence of the last section, 
explicit difference schemes which are high accuracy interpolation formulas seem 
most likely to succeed in obtaining accurate solutions of the transport equation. 
With this in mind, several new two-level formulas are now proposed and used to 
solve (2). These formulas give u,,,4; in terms of u at nodes on the time step s. 

I. Linear Interpolation Formulas. 


(7a) Ura = (1 ot P) Urs + PUr—-1,s (0 Ss Pp Ss 1) 
(7b) Ur si = (2 _ DP) Ur—1,s + (p TT 1) u,y-2, (1 < Pp > 2) 
(7c) Urq = (N+ 1 — P)Upnet (P—N)Urnie (NSEpSsn-+ il). 


These formulas are obtained in the following manner using Figure 5. If PP; is the 
characteristic through P so that up = up, and P, lies between B and C, then 
BP,:P,C = p:1 — p, and by using linear interpolation of u between B and C, 
formula (7a) is obtained. If the characteristic through P cuts the line s in P, 
between C and D, then CP,:P,.D = p — 1:2 — p, and linear interpolation of u 
between C and D gives formula (7b). A similar method may be used for values of p 
greater than 2. 

The general result for any value of p lying between integers n and n + 1, where 
n may be positive or negative, is given by (7c). The formulas are stable since the 
coefficients on the right-hand sides are positive and add to unity in corresponding 
pairs. 
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TABLE l(a) 
| Scheme 
} . 
Renee | — Carlson Central Difference 
r Values |— 
je} 4] 4 | 3 : + | 
i co. | 
rs | 12| 12 | 2 24 pm | 
| ic Sol 
~10 o | — |— |] — | oo; o | 
-9 wee paeepereks | —0.005 | —0.003 | 
—8 0 2 i ote Meer 0.007 0.002 | 
7 o | — | — | — _ | —0.028 | —0.024 | — 
—6 |S ee | — | — 0.029 0.015 | 0.00: 
—5 0 — |—]— —0.051 | —0.045 | —0.03: 
—4 0 | 0.001 | 0 — 0.012 0.006 | 
— | @ 0.03/0 | — | 0.053) 0.052 | 
—2 | 0 | 0.010 | 0 | — | —0.059 | —0.032 | —0.010 
—-1 | O | 0.017} 0 Seam. 0.006 | —0.024 
0 | O | 0.044/0 | 0.032} 0.008, 0.014] 0. 
1 | QO | 0.071 | 0.001 | 0.110) 0.166!) 0.166} 0.15: 
2 | 0 | 0.147| 0.011 | 0.241 | —0.105 —0.067 | —0. 
i oe 0.223 | 0.047 | 0.424 | —0.270 | —0.267 | —0. 
4 | 0 | 0.384 | 0.144 | 0.657} 0.021 | —0.022 | —0.02: 
5 0 | 0.545 | 0.338 | 0.934) 0.200) 0.213 | 
6 | 0.5 | 0.796 | 0.644 / 1.188 | 0.465 0.490 | 
7 | 1.0 | 1.046 | 1.044] 1.382] 0.842 0.816 | 
“Seca 1.311 | 1.488 | 1.498| 1.676 1.620 | 
9 | 2.0 | 1.575 | 1.906 | 1.532 | 2.403 2.388 | 
10 | 2.5 | 1.715 | 2.210] 1.487| 2.605 | 2.592 | 
11 | 3.0 | 1.856 | 2.323| 1.368] 2.631 2.652 | 
12 | 2.5 | 1.773 | 2.210| 1.216] 2.426 | 2.413 | 
13 | 2.0 | 1.691 | 1.906! 1.055 | 2.023 2.050 | 
14 | 1.5 | 1.428 | 1.488} 0.898 | 1.408 1.436 | 
15 | 1.0 | 1.166 | 1.044 | 0.753 | 0.817) 0.828 | 
16 | 0.5 | 0.873 | 0.644 | 0.625 | 0.424 0.440 | 
17 | O | 0.580 | 0.3388 | 0.514| 0.169 0.152 | 
18 0 | 0.385 |0.144| 0.419| 0.070 0.066 | 
19 0 0.190 | 0.047 | 0.339 0.019 0.012 
20 0 | 0.110|0.011| 0.274| 0.006 0.004 
21 0 0.031 | 0.001 | 0.219| 0 0 
22 0 0.015 | 0 | 017%; — — | 
2s. |. 0 0 0 | cima). — — | 
4% | O 0 0 | 0.110; — — | 
25 “-.. 0 0 0.087; — — 
wc ee @ 0 0 0.069; — — 
Bivch 0 0 | 0.054, — ~- 
28 0 — | — | 0.042; — — 
29 0 — — 0.033; — — 
30 0 — — | 0.006; — — 
31 o |} — | — | 002; — — 
Del | 7.298 | 2.907 | 12.308 | 3.000 2.723 
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TABLE 1(b) 
| Scheme 
| Friedrichs | Carlson Central 
, | Correct | Difference 
| Values _ 
ne? 1 ] 3 3 
s | 12 12 2 12 

0 —0.707 |  -—0.317 | —0.560 —0.477 | —0.671 

1 | 0.924 | -0.4389 | —-0.731 —0.388 |  —0.897 

2 | 1.000 | -0.495 | -—0.792 —0.269 | —0.996 

3 —0.924 | -0.475 | -—0.731 —0.128 | —0.938 

4 —0.707 ~0.383 | -0.560 | 0.022 —0.738 

5 —0.383 | 0.232 | —0.303 | 0.162 | —0.430 

6 | 0 —0.046 0 0.276 | —0.047 

7 | 0.383 0.146 0.303 0.346 | 0.350 

8 0.707 | 0.317 0.560 0.364 | 0.671 

9 0.924 | 0.439 0.732 0.327 | 0.897 

10 1.000 | 0.495 0.792 0.242 | 0.996 
11 0.924 | 0.475 0.732 0.121 0.938 
12 0.707 | 0.383 0.560 | -—0.016 | 0.738 
13 0.383 | 0.232 0.303 —0.150 0.430 
14 0 0.046 | 0 —0.261 0.047 
15 —0.383 —0.146 —0.303 —0.330 —0.350 
> lel 5.174 2.094 7.950 | 0.478 





II. Parabolic Interpolation Formulas. 
(8) Urex = —3(p — m)(n + 1 — p)trsine + (p —n + 1)(n +1 — p)tngs 
+3(p—n)(p—n+ 1) Upr-nue (nSpsnt 1). 


Referring again to Figure 5, if PP; is the characteristic through P, where P, lies 
between B and C, and if a parabolic interpolation formula incorporating the values 
of u at A, B andC is used to give u at points between B and C then u,,,4; is given 
by (8) with n = 0. If the characteristic through P cuts the line s at P, where P: 
lies between C and D, and a parabolic interpolation formula incorporating the values 
of u at B, C, and D is used to give u at points between C and D, then u,,.4; is given 
by (8) with n = 1, and so on for higher values of n. Finally, the stability of (8) is 
easily demonstrated by using methods described in Richtmyer [1], since the equa- 
tions are linear and have constant coefficients. Other stable parabolic interpolation 
schemes based on (8) are possible but they are unlikely to be more accurate than 
(8) with the original range of p stated. 
III. Cubic Interpolation Formulas. 


(9) Ureo = —3(p —n)(n +1 — p)(n +2 — p)trsins 
+3(n+2— p)(n+1— p)(pt+1—n)un, 
+ 3(p — n)(n + 2 — p)(pt+1 — n)u ns 
—t(p—n)(n+1—p)(ptil—n)uron,. (Spent). 
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From Figure 5, if P; lies between B and C, and a cubic interpolation formula based 
on the values of u at A, B, C, and D is used to give values of u between B and C, 
then u,,.4: is given by (9) with n = 0. If, however, PP? is the characteristic through 
P where P; lies between C and D, and a cubic interpolation formula based on the 
values of u at B, C, D and E is used to give values of u between C and D, then 
Ur,s+1 18 given by (9) with n = 1, and so on. Formula (9) is stable not only for the 
range of p stated but for the extended range n — 1 S p S n + 2, and so other 
cubic interpolation schemes based on (9) are possible. One other possible scheme is 
(9) together withn —-1 Spsin+2forn=--- —6, —3,0,3,6, --- . It is 
unlikely, however, that any of the other schemes will be as accurate as (9) with the 
original range of p stated and n any integer. 


4. Three-Level Formulas. So far the only three-level scheme discussed is the 
central difference formula. This gives u,,,4; in terms of u at nodes on the time steps 
s — 1 and s. Other three-level formulas suitable for limited ranges of values of p 
are now proposed. 
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Referring to Figure 6, PP; , GG, and HH, are the characteristics through P, G, 
and H respectively, thus up = up, , Uc = Uc, and Uq =Uz,. If P; lies between B 
and C, and a cubic interpolation formula incorporating the values of u at G, , B, 
H,, and C is used to give the values of u at points between B and C, then the value 
of u at P is given by 


_ (1 — 2p)(Q1 — p) 
l+p 


Urs = 





Ur,s—1 + 2(1 ns 2p )Ur.s 
(10) 
2p(1 — 2p) ' 
l+p 
In Figure 7, PP; , HH, , and IJ, are the characteristics through P, H, and I 
respectively, thus wp = up, , Uz = Ug, , and u; = uy, . If P; lies between B and C, 
and a cubic interpolation formula incorporating the values of u at B, H, , C and J, 
is used to give the values of u at points between B and C, then the value of u at P 
is given by 


+ 2pursa 21 ~~ 





r—l,s- 


2(2p — 1)(1 — 
won = — fee Ur.s + 2(1 ron Pp )Ur—1,2-1 


(11) 9 
+ 2(2p ant 1) ttt, a | oe Ur—2,s—1- 
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The stability of (10) and (11) for the range 0 < p S 1 can be demonstrated in the 
usual manner. 

Numerical calculations are now carried out using selected two- and three-level 
interpolation schemes to solve (2). The results are shown in Tables 2(a) and 
2 (b). The errors are shown in the last two rows where ).| ¢| is the sum of the 
moduli of the errors after a time 6 7 and >, | 1 | refers to the errors at a later stage 
in the computation when the profile has been transported over a further time 
6 + . In the case of the three-level formula (11), after a time 36 + the sums of 
the moduli of the errors are still only 0.660 and 0.026 for the ‘‘roof top” and “‘sine”’ 
curves respectively. The results shown in Tables 2(a) and 2(b) are for values of p 
lying between 0 and 1, but in the case of the two-level schemes they may be inter- 
preted for values of p outside this range. For example, the figures for p = 3 refer 
also to p = n + }$ if the profile is moved on a further 12n intervals of x. 


5. Interpolation Formulas and Finite Difference Schemes. In view of the form 
of the transport equation, a close link might be expected between interpolation 


1 


on 
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TABLE 2(a) 





Scheme 





| Parabolic 


Cubic | Three-Level | Three-Level 
Correct | Formula (8) | Formula (9) | Formula (10) | Formula (11) 














: Values 

| “SC WW WING oc 
s ict 12 | 24 

= 
4 | 0 0 0 0 0 
—3 | 0 0.001 0 0 0 
—2 | 0 —0.004 0 0 0 
—) | 0 —0.004 0 | 0 0 
Oy P< 0.020 | 0 0 0 
— 0 0.024 0.003 0 | 0 
4 0 —0.045 —0.005 0 0 
as | Ss ~0.126 | -0.032 | -0.003 | 0 
ie 0 —0.068 —0.026 | -0.027 | -0 
5 0 0.165 0.121 0.063 | 0. 
6 0.5 | 0.502 | 0.471 0.471 | 0. 
7 ie | 0.961 | 0.961 1.003 | 1. 
8 15 = | 1.591 | 1.505 | 1.498 iy 
9 2.0 | 2.248 | 2.068 2.007 hy 
10 2.5 | 2.651 | 2.554 2.555 2. 
il 3.0 | 2.682 | 2.757 2.872 2. 
12 2.5 2.433 | 2.554 2.557 | 2. 
13 20 | 2.006 | 2.068 1.992 | 1. 
m3 >» £8 1.452 1.505 1.501 | 1. 
15 1.0 0.876 0.961 0.996 | 1. 
16 0.5 0.421 0.471 0.472 | 0. 
17 0 0.156 0.121 0.063 0. 
18 0 | 0.043 —0.026 —0.028 —0. 
19 0 | 0.008 —0.032 0.003 0 
20 0 | 0.001 —0.005 0 0 
21 0 0 0.003 0 0 
22 0 0 0 0 0 
23 0 0 0 0 0 
24 0 0 0 0 0 
25 0 0 0 0 0 
6 0 0 0 0 0 
> lel | 1.838 1.007 0.509 0. 
> | al | 2.644 1.169 | 0.660 0 





| 











oo 


28 


—_ 


287 


.432 





formulas and difference schemes used to solve (2). This is best illustrated by means 
of an example. Consider the problem of evolving a finite difference replacement of 
(2) which makes use of the points P, B, C, and D in Figure 5. Taylor expansions 


about the point B give 


(12) 


(13) 


ot 


Ur.et+l = Ur,s + At () 


Us-it,s = Urs — Az ( 


au 
0 


2 au 
1 —- 
*) + , (Az) (=). 
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TABLE 2(b) 
Scheme 
Parabolic Cubic Three-Level | Three-Level 
3 Correct Formula (8) | Formula (9) | Formula (10) Formula (11) 
Values 
p | 3 } t : 
s | 12 12 24 8 
0 —0.707 —0.733 —0.702 —0.706 —0.706 
1 —0.924 —0.933 —0.917 —0.922 —0.923 
2 —1.000 —0.992 —0.993 —0.$39 —0.999 
3 —0.924 —0.900 —0.917 —0.923 —0.923 
4 —0.707 —0.670 —0.702 —0.706 —0.706 
5 —0.383 —0.338 —0.380 —0.382 —0.382 
6 0 0.044 0 0 0 
7 0.383 0.420 0.380 0.382 0.382 
8 0.707 0.733 0.702 | 0.7 0.706 
9 0.924 0.933 0.917 | 0.922 0.923 
10 1.000 0.992 0.993 0.999 0.999 
11 0.924 0.900 0.917 | 0.923 0.923 
12 0.707 0.670 0.702 | 0.706 0.706 
13 0.383 0.338 0.380 | 0.382 0.382 
14 0 —0.044 0 0 0 
15 —0.383 —0.420 —0.380 | —0.382 —0.382 
> lel 0.460 0.074 | 0.016 0.014 
> | a | | 0.914 0.144 | 0.026 0.014 
and 
ou 2 au 
(14) Ur—2,3 = Ur,s — 2Ax + 2( Az) =a.¢ . 
ax re Ox? r.8 


ou 
The value of (%’) 


r,s 


is obtained from (12) and (%*) 


r,8 


by climinating( 2) 
dx" 


r,s 


from (13) and (14). These values are then substituted into (2) to give 


(15) 


Ure = (1 — *?) Ur,e + 2pUr-1,s — P us 8 


2 


The truncation error in (15) is dominated by the term }( At) (F *), » neglected in 
(12) and since by differentiating (2) the result 


(16) 


is obtained, it follows that the principal part of the truncation error is }p" (Ar) 


ot 


72 Ju 
Ox? 


ect 


This is the standard finite difference approach which can, however, be improv ef in 
the following manner. Replace equation (12) by 


(17) 


tats = tra +t (SH) + 3 (At)’ ($).. ’ 





ie uw é -~ 
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which, on using (2) and (16) becomes 


(18) trot ty pas (% =) +4 pi(acy (Zt) 


If () and (=) are now eliminated from (13), (14), and (18), the parabolic 
interpolation formula (8) with nm = 1 is obtained with truncation error 


ip(1 — p)(2 — p)(A2)* 2 


This is a distinct improvement over the previous finite difference formula (15), 
and in particular if p is close to unity, the interpolation formula is expected to be 
specially accurate when used to solve (2). If p = 1, of course, the theoretical solu- 
tion of the interpolation formula (8) with n = 1 is the same as the theoretical solu- 
tion of (2). However, as it is intended to use the results of the present investigation 
to solve the general transport equation (1), the exact correspondence of the theo- 
retical solutions of (2) when p = 1 can really be ignored. This example illustrates 
the fact that the best finite difference formula for a given set of points used to solve 
(2) is an interpolation formula. This is because each derivative with respect to a 














TABLE 3 
Scheme | ctor ony | Truncation Error 
| | a3 
Friedrichs | (4 MA» ‘(a2 - 
CarlsonOSpZi1 | _ (a) 4p(1 — p)(ax)* 2 
Carlson p > 1 (5b) 3p(p + 1)(4z)* fe 
Central Difference | (6) —ip(1 — neal 
Linear Interpolation (7c) —3(n — p)(n — p+ 1)(Ax)’ — we > 
3 
Parabolic Interpolation (8) |—3(n — p — 1)(n — p)(n — p + 1)(Az)’ = 
Cubic Interpolation (9) |A(n —p— 1) (n— p) (n-—p+1) 
yy 4 a u 
| n — p + 2)(Az) ari 
—4p'(1 — p)(1 — 2p)(az") 2% 
Three-Level I (10) 3 : 
4 S u 


Three-Level II (11) | tp(1 — p)*(1 — 2p)(Ax) = 
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co-ordinate is a constant multiple of the corresponding derivative with respect to 
the other co-ordinate, thus the Taylor expansions can all be expressed in terms of a 
single variable. Elimination of the maximum possible number of derivatives with 
respect to this variable leads to an interpolation formula. 


6. Truncation Errors. For purposes of comparison, the truncation errors asso- 
ciated with the finite difference schemes considered for solving (2) are given in 


Table 3. The errors quoted are Az times the errors as defined by Richtmyer 
{1, p. 19]. 


7. Conclusions. The calculations carried out in the present paper, using existing 
stable finite difference schemes in turn to solve the simplified transport equation 
(2), vary considerably in accuracy. The central difference formula (6) is most 
accurate with Carlson’s scheme for |p| < 1 next in order of merit. Carlson’s 
implicit scheme for |p| > 1 is very poor, particularly for large values of | p|. 
This is illustrated in Figure 8 where the part of the truncation error depending on p 
is shown as a function (2) of p. It can be seen that the maximum value of the trun- 
cation error when 0 S p & 1 is one-eighth of the minimum value when p > 1. In 
fact, the authors believe that implicit schemes can be abandoned as a means of 
obtaining accurate solutions of the transport equation. 

New explicit schemes, derived as interpolation formulas, are next used to solve 
(2) and a considerable improvement in accuracy is obtained, particularly for 
schemes such as (9), (10), and (11), which are cubic interpolation formulas with a 
very small truncation error. The error in any numerical solution of (2) takes the 
form of a smoothing out of the initial profile together with, in most cases, a super- 
posed stable oscillation. 
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It cannot be emphasized too strongly, however, that schemes which successfully 
solve (2) do not necessarily give comparable accuracy when used to solve (1), 
where V is a function of z, ¢, and u. On the other hand, difference schemes which fail 
to give accurate solutions of (2), can hardly be expected to be more successful when 
used to solve (1). The main difficulty in solving (1) numerically arises from the 
fact that the characteristics are curved and the distance BP, (Figures 1, 2, 3, 5, 6,7) 
is no longer given simply by V At or pAz. It must be found by integrating the equa- 
tion 


(19) & — V(z,t,u) = 0. 

If, in the case of curved characteristics, BP, is now expressed as p’ Ar, any one of 
the interpolation formulas proposed in the present paper may be applied directly 
with p’ substituted for p. The value of p’ is, of course, in general different at each 
node. 

In deciding the values of Az and At for a given calculation, Az is first chosen to 
represent adequately the initial profile. The time step Af is then chosen so that 
BP, is large enough for the calculation to proceed without too many interpolations 
but not so large that the positions of P,; , obtained from (19), are too much in error. 
We hope to examine in detail at a later date the general problem of integrating 


(1). 


8. Acknowledgements. The calculations were performed on the Deuce computer 
of Glasgow University, and the authors are indebted to Dr. D. C. Gillies, Director 
of the Computing Laboratory, for permission to use the machine and to Miss A. 
Low for help in running the programs. 


Mathematics Department 
Queen’s College 
Dundee, Scotland and 


Mathematics Department 
St. Salvatcr’s College 
St. Andrews, Scotland 


1. R. D. Ricutmyer, Difference Methods for Initial-Value Problems, Interscience Publishers 
Inc., New York, 1957. 

2. K. O. Friepricus, ‘‘“Symmetric hyperbolic linear differential equations,’’ Comm. Pure 
Appl. Math., v. 7, 1954, p. 345-392. 

3. B. G. Cartson, Unpublished Los Alamos Report, 1953. 

4. H. B. Kevier « B. Wenprorr, “On the formulation and analysis of numerical methods 
for time-dependent transport equations,’’ Comm. Pure Appl. Math., v. 10, 1957, p. 567-582. 

5. J. S. Matxus & G. Witt, ‘“‘The atmosphere and the sea in motion,’’ Rossby Memorial 
Volume, Rockefeller Institute Press, New York, 1959, p. 425. 








Quadrature Formulas for Infinite Integrals 


By W. M. Harper 


1. Introduction. Since the advent of high-speed computers, ‘mechanical’ 
quadratures of the type 


(1) [ wane) de ~ YH f(a) 


have become increasingly important. The only quadrature generally available for 
the case b = —a = ~ is the Hermite-Gauss formula although the Laguerre- 
Gauss formula can also be used if f(z) is an even function of x. The latter would, 
however, require computation of twice the number of ordinates for a corresponding 
degree of precision and would therefore rarely be preferred. In either case the in- 
tegrand is supposed to behave like the product of an exponential function and a 
polynomial. For purely algebraic integrands it would appear to be more appropriate 
to use a quadrature based on an algebraic weight function even though the degree 
of the polynomial approximation to f(z) is limited. 

In this paper, formulas of type (1) are derived with weight function w(x) = 
(1 + 2)“ for the range b = —a = ~. In a modified form they are shown to 
be superior to the Hermite-Gauss and Laguerre-Gauss quadratures for a particular 
class of statistical integrals. 


2. Derivation of Quadratures. In the quadrature formula 


(2) f (1+ 2)'Y(2) de = > H; f(a;) + En, 


the abscissas a; will be the zeros of the nth degree polynomial ¢,,.(2) which satis- 
fies the orthogonality condition 


(3) / (1 + 2?) bm (2) bn e( x) dx = 0, (m #¥n,m+n < 2k+1). 


By standard methods given for example in [2], [4], it is easily shown from (3) 
that the orthogonal system of polynomials is given by the Rodrigues formula 





— a (2k — 2n + 2) 2\k+1 a” 2\n—k—-1 


(n<k+1) 
where the standardizing constant is chosen to make the coefficient of x” unity. 


By direct manipulations with (4) and repeated use of Leibnitz’ formula, the re- 
currence relations (5)—(10) are easily established. They are 


n(2k — n + 2) 
k — 2n + 1)(2k — 2n + 3) on-1,4(2), 


(6) (1+ 2°)onatz) = (2k — m+ 1)xbaa(x) — (2k — In + Mobasia(z), 


(5) On41,4(2) _ hn e( 2) ar (2 
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n(2k — n + 2) 


2 ’ on 
(7) (1 + 2)onu(z) = nepnalz) + Ss 





On—i4(2), 
2k —2n+3 
(2k — n + 2)(2k — n +3) 
“[{(4k — 2n + 3) + (2k — 2n + 1)2"} bna(z) 
— (2k — 2n + 1)(1 + 2° )onea(z)], 
21+ 2 )onu(x) = [na* — (2k — n + 2)onn(z) 
(2k — n + 2)(2k — n+ 3) 
2k —2n+3 
(10) Tn(Z) = (2k — n+ Lbau(z) — (2k — 2n + I) bnna(z). 


The polynomial system can now be extended to include values of n excluded in (4). 
For n > k + $ however, complex zeros make their appearance so that no useful 
quadratures are available for this range of n. 


It is similarly easily shown that ¢,.(z) is a solution of the differential equation 





On e4a( 2) bal 
(8) 


(9) 
bp 





on aui(z), 


(11) (1 + 2°)y” — 2kay’ + n(2k —n+ 1)y =0 
whence the relation 

pen t fl ,r(k — n + $) y —k-1/2/ 
(12) Gn e(X) (5) "—TeS C. (iz) 


can be established where in the notation of [1], C,(z) (designated by P,(z) in 
[7]) is the Gegenbauer or ultraspherical polynomial of degree n and parameter \. 
Relations with Legendre functions can also be established, namely: 
na(z) = (—1)"" x"? lim 
sk 
é [ren I(s ~™ @ + #) 
y T(2s — n + 2) 
where P,”(z), in the notation of [1], is the associated Legendre function of the first 
kind with parameters » and v, and 


(14) dna(z) = AP (k — m+ HL + YO PTA + 2] 


where P,”(z) is the associated Legendre function of the first kind with definition 
suitable for the cut in the real axis from z = —1 toz = 1. The limit in (13) caters 
to integer values of k (see [3]). 

The weight coefficients and error term in (2) can be determined by standard 
methods with the results 


(13) 








cosec sx(1 + 2ynnpth is) | 








‘ 3)]? ’ a 
(15) H; = 2*°""n! Ren tH (1 + aj) ‘lonn(a;)I”, 
Od ee ko 


gent nir(k —n + $)P (2n) ( 
De | 
Gk — In + DOn)ITOk— na! n<k+! 
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The restriction on n is necessary to ensure convergence of the error estimate but 
does not ensure that a close upper bound to the actual error can be obtained (see, 
for example, [2]). 

For practical purposes a more convenient form of the quadrature is 


(17) [_ He) az ~ & Ks fla); 
c) j= 
here the weight coefficients are given by 


(18) K; = H,(1 + a?)*™. 


The values of a; and K; for four- and six-point formulas for some integral values 
of k are given in Table 1. 

The right-hand side of (17) is a function of k as well as of n; for a given value 
of n, therefore, there will be a value or values of k depending on f(x) which will 
give the “best” approximation to the integral on the left. The determination of 
such values and the corresponding parameters appears to be too formidable a 
task for practical applications. For the special cases k = n — 1, k = n, however, 
solution of (11) with z = cot @ enables ¢,.(z) to be obtained in the forms 


(19) Gn,n-1(2) = cosec” (are cot x) cos (n are cot z), 
(20) On.n(x) = (n + 1) cosec”™ (are cot x) sin [(m + 1) are cot z]. 


The zeros are now simple cotangents and the weight coefficients H ; assume simple 
trigonometric form; the resulting quadratures can be written as 


(21) | are e (k=n-—1), 


j=l 
ze), (k =n). 


These formulas can also be deduced from the Chebyshev-Gauss quadratures 
' ~ _ (25 = 1) 
(23) [ (1 — #)*o(u) dy ~ZD gfe eae 


(24) [ (1-7) gy) dy ~ — 7 de (co 8 7.) 


by the substitutions y = x(1 + 2°)”, g(y) = f(z). 








(22) i (1 + a’) f(a) dx ~ are ie 2, § (cot - 


3. Practical Application. An example of a useful application for the quadratures 
is the evaluation of integrals arising in the determination of the statistical distribu- 
tion of the ratio of two quadratic forms in normal variates. If the quadratic forms 
are independent mean half-square successive differences based on sample sizes of 
p and q respectively, one of the integrals which require evaluation can be written 
in the form 


c p—l q—1 
‘ I(z) -_ | (1 af. a)" T] (a, + a)? Il (1 + bee” 4p er dx, 
(25) 0 r=2 s=l 


(p even), 





TABLE 1 
Abscissas and Weights for Quadrature (17) 











— 























A.n=4 
k +a; K; 
3 0.41421 35624 0.92015 11845 
2.41421 35624 5.36303 41227 
4 0.32491 96962 0.69465 18830 
1.37638 19205 1.81862 22399 
5 0.27618 30252 0.58086 65620 
1.06005 79874 1.17945 11502 
6 0.24436 83118 0.50932 47880 
0.89298 76737 0.90816 46087 
7 0.22150 78137 0.45903 94023 
0.78587 59159 0.75578 97944 
| 
8 | 0.20405 97869 0.42121 27662 
| 0.70979 86678 0.65698 70999 
| 
9 | 0.19017 76238 0.39142 46836 
0.65220 46710 | 0.58705 73261 
| 
10 | 0.17879 14705 0.36717 90805 
0.60665 77372 0.53455 96626 
B. n = 6 
k +a; K; 
5 0.26794 91924 0.56119 14763 
| 1.00000 00000 1.04719 75512 
| 3.73205 08076 7.81638 89333 
6 0.22824 34744 0.47217 91694 
0.79747 33889 0.73421 88392 
| 2.07652 13966 | 2.38399 35955 
| 
7 0.20219 30919 | 0.41550 76425 
0.68370 47228 0.58969 00381 
| 1.57850 04858 1.44716 80133 
8 | 0.18342 80037 0.37535 93234 
0.60816 30047 0.50404 67421 
| 1.31884 38384 1.06492 43997 
i 0.16907 35256 | 0.34499 40643 
0.55326 32106 | 0.44635 57833 
1.15411 46518 0.85743 60559 
10 0.15763 63749 0.32098 68394 
0.51101 94490 0.40432 69556 
0.72680 65190 


1.03809 74230 
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TABLE 2 
Comparison of Quadrature Formulas in Evaluating I(i) 

















Quadrature No. Abscissas Result E X 10° 
Series — 1.2106 5423 — 
Algebraic, k = 5 6 1.2106 4384 1039 
Algebraic, k = 6 6 1.2106 5381 42 
Algebraic, k = 7 6 1.2106 5415 8 
Algebraic, k = 8 6 | 1.2081 0423 25 5000 
Algebraic, k = 9 6 1.2025 0816 81 4607 
Algebraic, k = 10 6 1.1942 4044 164 1379 
Hermite 6 | 1.1610 8623 495 6800 
Hermite 8 | 1.1879 0738 227 4685 
Hermite 10 =| 1.1994 3337 | 112 2086 
Laguerre 6 iB 


.1674 2007 | 432 3416 





where the a, and b, are constants. In order to compare methods (25) was evaluated 
by various quadratures for the case p = 4, g = 3, z = 1 when the test integral 
becomes 


ra) = [ a+ey|(bv2+2) @v2-242) 
, {5 (7 — 24/2) + a {5 (13 — 24/2) + zy" dr. 


The quadrature (17) was applied for the values k = 5(1)10 using six abscissas 
in each case. The Hermite-Gauss quadrature was used with six, eight and ten 
abscissas, and the Laguerre-Gauss formula for six abscissas (which requires the 
same number of evaluations of the integrand as the other formulas for twelve 
abscissas but which is of degree of precision eleven as against twenty-three for the 
others). The abscissas and weights for the Hermite formula were taken from the 
values tabulated in [6] and those for the Laguerre method from [5]. The results 
together with the correct value of J(1) determined by a series method are tabu- 
lated to eight decimal places in Table 2 which also shows the errors of the methods. 

The table shows the superiority of the “algebraic” quadratures over the Her- 
mite and Laguerre formulas for this integral; even the use of ten abscissas for the 
Hermite quadrature leaves an error much greater than the algebraic quadratures 
with only six abscissas except for the case k = 10. The best algebraic quadrature is 
for k = 7 but the advantage over those for k = 5 and k = 6 is too small to com- 
pensate for the simplicity of the latter two cases when used in the equivalent forms 
shown in (21) and (22) respectively. In addition, the quadrature (22) evaluates 
I(1) correctly to eight decimal places for n = 8 as does (21) for n = 9. 


(26) 
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A Modification of the Runge-Kutta 
Fourth-Order Method 
By E. K. Blum 


1. Introduction. Consider the system of n first-order ordinary differential 
equations, 


(1.1) yx’ = felt, w(t), «~~ , ya(2)), hie 3) 00 
with the initial values, 
(1.2) yi.( to) = a. 


Under suitable conditions on the f;, a unique solution of (1.1) satisfying (1.2) 
exists for some interval, 4; < ¢ < b. For example, it is sufficient that the f; be con- 
tinuous and satisfy a Lipschitz condition in some neighborhood of the initial point, 
(to, @1,*** , @,). We shall assume that such conditions obtain, so that the initial 
value problem (1.1), (1.2) has a unique solution. 

To simplify the notation, we define yo = ¢ and fy = 1. We now let y be the 
vector, (Yo, ¥i1,°**» Yn), and f the vector-valued function, (fo, fi,---, fn). 
The initial value problem can then be written as 


(1.3) y’ = f(y), 
(1.4) y(t) = a. 


The Runge-Kutta fourth-order method for the numerical solution of (1.3), 
(1.4) yields approximate values, y;, of y on a finite set of points, t; = .)» + jh, 
j = 1, 2,---, m. It is usually summarized in formulas (1.5)-(1.9) below, which 
specify the calculations to be carried out for each integration step; i.e. for each 
value of j. 


(1.5) ky = hf(ys), 

(1.6) ke = hf(ys + k/2), 

(1.7) ks = hf(y; + ke/2), 

(1.8) ky = hf(y; + ks), 

(1.9) Yiar = Ys + (hy + Bho + By + ky) /6. 


A variant of this method was derived by 8. Gill [1]. The two advantages of Gill’s 
variant are (1) in automatic computers, it requires 3n + B storage registers 
whereas the Runge-Kutta formulas as given above, require 4n + B, where B is 
some constant; (2) the computation can be arranged so that rounding errors are 
reduced appreciably. In the present paper, we shall show how, by means of a fairly 
simple modification of (1.5)—(1.9), both of these advantages can be made to accrue 
to the classical Runge-Kutta method. All the constants in this modification are 
rational, whereas Gill’s variant contains some irrational constants. The modifica- 
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tion is achieved by extracting from Gill’s method its main virtue, thé rather in- 


genious device for reducing the rounding error, and applying it to a rearrangement 
of (1.5)-(1.9). 


2. The Exact Modification. In an automatic digital computer, real numbers are 
replaced by what von Neumann and Goldstine [2] call ‘digital numbers,” that is, 
by real numbers rounded to a prescribed number of digits. Further, exact arith- 
metic operations are replaced by “pseudo-operations” since results must be 
rounded. The main advantage of the modified Runge-Kutta formulas to be pre- 
sented in Section 3 is that they reduce considerably the rounding error arising 
from the unavoidable use of digital numbers and pseudo-operations. The saving of 
n storage registers is a secondary consideration in large computers. The same is 
true of the Gill variant. 

In this section we shall present a preliminary version of the proposed method. 
We shall refer to it as the “exact modification” since all operations will be assumed 
to be exact operations on real numbers. The form of the exact modification will 
demonstrate clearly how the saving of n storage registers is effected. 


Using vector notation, as in (1.3)—(1.9), we can write the exact modification in a 
recursive form as follows: 


( 2 = Yi,» 
(2.1) es = Yj, 
Po = hf(z), 
( 2, = zo + Py/2, 
(2.2) is = Pe, 
P, = hf(a), 
2 = 4 + P,/2 — @,/2, 
(2.3) gq = 9/6, 
P, = hf (z2) = P,/2, 
23 = 2 + P2 ’ 
(2.4) 3 = q — P2, 
P; = hf(zs) + 2P2, 
(2.5) Yjun =a = 2% + qs + P;/6. 


(Strictly speaking, each of the vectors, z; , gi: , P: , should have a second subscript, 
j, to indicate that the sequence (2.1)—(2.5) is repeated for each step of the solu- 
tion. This subscript has been dropped for reasons of economy, just as the subscript 
which indicates the components of the vectors has been dropped.) 

THEOREM 1. The exact modification, (2.1)~(2.5), is equivalent to the classical 
Runge-Kutta method and requires only 3n + B storage registers. 

Proof. To show that (2.1)—(2.5) is equivalent to (1.5)—(1.9), we first observe 
that Py) = k,. Then z = y; + k:/2, which implies P; = k.. Since q = hk, 
it follows that z. = (y; + k:/2) + ke/2 — k,/2 = yj + keo/2. Thus, 


P, > ks = k2/2 
and gq: = k,/6. From (2.4), it now follows that z; = (y; + k2/2) + (ks — ke/2) = 
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Y3 os ks ; and qs = k;/6 - (ks = ke/2), whence P; at ks _ 2Qkz _ ke ‘ Combining 
these expressions in (2.5), we get 


ke 


(uy + ha) + (B+) +2 (he + ke — be) 


Yin 


Sosa =e + ; (Jey + he + 2ha + ke). 


The order of computation of the components of the vectors 2:41, @i4i1, and 
Pii1, i = 0, 1, 2, 3, should be as follows. First, compute the components of 2:4; 
and q;4: together. Each such component involves only the corresponding com- 
ponent of z; , q;, and P;. Hence, as each component of z;4; and g;4: is computed, 
it can be placed in the storage occupied by the corresponding component of z; 
and gq; , respectively. After all components of z;4; and gi;; have been computed, 
the components of P;,; can be computed, replacing the corresponding components 
of P; in storage. Thus, 3n + B storages suffice. 


3. The Finite-Precision Modification. In this section we shall consider the 
rounding errors which arise in actual computation when digital numbers and pseudo- 
operations are used in (2.1)—(2.5). We shall adopt the notation of [2] for digital 
numbers, denoting a digital number by a letter with a bar over it, and similarly 
for vectors; e.g. 7; is a vector having digital numbers as components. However, 
we shall not introduce special symbols for pseudo-operations. Instead, we prefer to 
write all formulas with exact operations and introduce special terms to denote the 
rounding error caused by the pseudo-operations. Besides the usual arithmetic 
operations, we require two “shifting” operations. These are best described in- 
formally. 

For the remainder of this section, let us assume that a digital number is repre- 
sented by a sequence of s decimal digits, and that the decimal point is at the ex- 
treme left. (The first digit immediately to the right of the decimal point is said 
to be in position 1.) For m a non-negative integer, we define the operator, R,, , 
(“shift right m places”) as follows. If 7 is a digital number, then R,.g is the digital 
number obtained by shifting the digits of 7 m positions to the right, “rounding off’’ 
the digits shifted into positions s + 1, --- , s + m, and inserting zeros into posi- 
tions 1, --- , m. The usual method of rounding off is to add (if 7 = 0) or subtract 
(if 7 < 0) the digit “5” in position s + 1, and then drop the digits beyond position 
s. The operator, L,, , (“‘shift left m places”) is defined similarly. Thus, Lj is the 
digital number obtained by shifting the digits of 7 to the left m places, dropping 
those digits which are then to the left of the decimal point, and inserting zeros into 
positions s — m + 1, --- , s. When applied to vectors, R,, and L» are considered 
to operate on each of the component digital numbers. If y is a real number, then 
we define Ray = 10°"y and Lay = 10”y. 

To briefly motivate the formulas to be given below, let us consider the exact 
modification, (2.1)—(2.5). If this procedure were carried out with digital numbers 
and pseudo-operations, an analysis of the rounding error would show that, under 
suitable conditions on the partial derivatives, df,/dy; , the main source of error is 
in the computation of the z;. The error there arises from the fact that g; and P; 
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are usually smaller than z; in magnitude, since they are of order h. This means that 
either g; and P; must be computed with fewer significant digits than z; , or else a 
shift right must be performed before g; and P; can be added to z; . In either case 
an appreciable rounding error is incurred. The procedure explained below reduces 
this particular error. We shall refer to it as the “finite-precision modification” of 
the Runge-Kutta method to emphasize that it is designed for actual computation 
with digital numbers having a “finite precision” of s places. The finite-precision 
modification is derived from the exact modification by introducing the special 
quantities, r; , as in Gill’s formulas. To facilitate the error analysis, we shall write 
the finite-precision modification first with real numbers and exact operations, 
(3.1)-(3.5), and then with digital numbers and error terms, (3.1)—(3.5). 


2 = Vi» 
(3.1) Lingo = Lunds; ; 
LmPo = (Lh) f(z), 
( 11 = Ru($LmPo — Lago), 
ve =zat+n, 
Lng ae 3L al <« (4L Po r Lino) ; 
LP; = (Lwh)f(a), 


2 = Ra(4(LaP = Lm) ); 
: = 2; + T2, 
(3.3) Ve Le = -Ln — jhe + it. 


LmP2 = (Lmh)f(z2) — iLmPs , 


_ R»(LnP2), 
=z + T3, 
(3.4) aa ae —Lats + Linge ’ 
LP; = (Lwh)f(zs) + 2LnP2 , 


ts = Ra($LaPs + Langs), 
(3.5) | Yu=wa=—“2+n%, 
Lin4,541 3 [Lats - (4LnPs + Lnqs)). 


Remark 1. Regarding all operations in (3.1)—(3.5) as exact, we can replace 
R,, by 10°” and L,, by 10”. A straightforward calculation then shows that (3.1)- 
(3.5) is equivalent to (2.1)—(2.5). 

Remark 2. The quantities, r;, are redundant if all operations in (3.1)—(3.5) 
are considered to be exact. They play a significant role only when digital numbers 
and pseudo-operations are introduced. 

Remark 3. The r; require only one additional storage register rather than n. 
The order of computation should be as follows. First, compute together the com- 
ponents of r;, z;, and Lmg;, as indicated by the inner brackets. Then the com- 
ponents of L,.P; can be computed. Since the computation of a component of 2; 
and Lng: requires only the corresponding component of r; , and since r; is not used 
after the ith stage, one storage suffices for all components of all r; . 

Remark 4. It is obvious that the quantity, Lngqs; , is always zero if exact opera- 
tions are used. For pseudo-operations this is not the case. However, Lng can 
always be taken to be zero to start the computation. 
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In practice, the finite-precision procedure, (3.1)—(3.5), would be carried out 
with digital numbers and pseudo-operations. The operation “+” would be ex- 
ecuted as pseudo-addition, and both R,, and L,, wouid be executed as shift opera- 
tions rather than as multiplications. To analyze the rounding error, it is convenient 
to rewrite (3.1)—(3.5) in a mixed form, (3.1)—(3.5), involving digital numbers, 
exact operations, and error terms. It is to be understood that the digital numbers 
are thereby treated as real decimal numbers having zeros in all positions beyond 
position s. The effect of pseudo-operations is shown by the presence of a single 
error term denoted by an expression of the form, e(u). 











( a=9; 
(3.1) ie = Lau, 

LimPo — (Lmh) f (Zo) + e(LmPo), 

7, = Ra(4Lngo — LnPo + e(r:), 
(3.2) 4A=a+in, 
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5 
cy 

I 


3Lmht — ($LmPo — Lmgo) + e(Lmg), 
= (Lwh)f(%) + e(LmPs), 

= Rn(}(LmPi — Dmqi)) + e(r2), 
=Aa+h:, 

—LmFz — 4Lmqi + $LmP, + e(Lmqe), 
= (Lwh)f(%) — 4LmP: + e(LaPs), 


i, = Rn( LmP2) + e(rs), 
B=-hTh, 
LmQs = —Lnits aa Linge pub ls 
(Lmh)f(%3) + 2LmP2 + e(LmP3), 
Fy = Ru(§LmPs + Lngs) + e(1s), 
(3.5) Yj — 24 nian 23 + Ms, 
Tnti.iss = 3{Lnfa — (4LwPs + Dogs)] + e(Lnge)- 
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Remark 5. It is seen that e(z;) = 0 for i = 0, 1, --- , 4 because the pseudo- 
operation of addition gives the same result as the exact operation. This is true 
because of our assumption that in all digital numbers the decimal point is in a 
fixed position. If “floating-point” numbers are used, the pseudo-operation of addi- 
tion can introduce a rounding error. We shall discuss this in the next section. 

We are now in a position to estimate the rounding error in (3.1)—(3.5). After 
some preliminaries, we shall formulate the results as Theorem 2 and its corollary. 

For any quantity, uw, we define «(u) = &@ — u; ie. e(u) is the total rounding 
error in u. We are interested in e(y;). However, it will turn out that the quantity, 

Rn 


i; = 9j i gz Lm M5 





is a better approximation to y; than is 9; . Thus, we shall consider e(y,;*) instead, 
where 


* RK. 1 
ym = — = bm Oi eh ~ 3 - 
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By remark 4, q4; = 0, so that y;* = y; . Hence, 


yi = 97 — ey#*). 
We note that 


e(y;*) = €(y;) — Am (Lm Qu) 


It is convenient to deal with the norm of a vector, u, which we define as 


|| a || = max | us|, 


where u; are the components of u. For a matrix, A, with elements, ay , we define 
|| A || = max {> | au |}. 
7 7 


In particular, we shall be concerned with matrices for which ay = df;/dy% , 
where the partial derivatives are evaluated at different points for each 7 and k. 
A matrix of this type will be denoted by the symbol, “J.” 

THEOREM 2. For any of the quantities, u, computed in (3.1)—(3.5), let the error 
term, e(u), be subject to the condition, 


(i) lle(u) | < 107. 


Let the bounds on the partial derivatives, df;/dy, , be such that for any matrix, J, 
(ii) | J || = L. 
Let hh = 10°", 0 < m < 8. Then the total rounding error in y;* incurred in one in- 
tegration step is not greater than 2M10*™ in absolute value. 
Proof. From (3.2) and (3.2) we obtain 
e( 2) = €( Zq) + Ra(4e(LmPo) = €( Limo) ) + e(ri). 
From (3.1), (3.1), if we assume that h = h, we have 
e(LmPo) = Lm(h)(f(9i) — f(yi)) + e(LmPo). 
Now, for each component, f; ,7 = 0,1, --- ,, we have 
" ~. Of; 
§(Gs) — fus) = DL egy), 
k=O OY jk 
or, in matrix notation, 
f(Gi) — flys) = Se(ys). 
This gives 
(Lm Po) = Lm(h) Je yj) + e( Lm Po), 
(3.6) h ae 
e(z1) = (yj) + 5 Jelys) — Ru Lm Gj) + > ChLe Po) + e(r:). 
Proceeding in this way, we obtain 


h 
2 


_— 


(3.7) e(z2) = ( +-J+ : rau) + e(r2) — =e(ri) +2 


~ 


t 
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where 
2m AS ely) + AP e(Lm Pr) — Se Lm gr) — ME" Sle Lm aus) — eC Lm Pa) 
le? os 1 
(zs) = {1 +hJ + 5 J+ a7 ys) — 5 ers) + e(re) + e(rs) 
(3.8) 
+ Rm ¢(Lm Ps) — 32 e(Lm qu) + hers) — 5 Je(ni) + he, 
| Rn 
(3.9) e(ys41) = e(ys) — €( Lm Gas) + hJv + .y W + e(r), 
where 
v = $(e(ys) + 2e(z) + 2e(z2) + €(2s)), 
and 


W = e(LmPo) + 2e(LmP1) + 2e(LmP2) + e(LmP3) + 6e(Lmge) — 2Rme( Dmg). 
Using (3.6)—(3.8), we obtain 


h h h’ Rn 
v= (1 + 57 + git rs) €(y5) a 3 €( Lim 443) 
(3.10) . ; 
+ 75 ot) + 5 etre) + 5 ers) + wi + mw, 
where 


~~ a (e( Lm Po) + €(Lm P1) + e(Lm P2)) — = er Sy, 


w= 2 +E I(e(m) — Sens) + 2). 


From (3.5) and the fact that q,; = 0, we obtain 
(3.11) €(LimQs,j41) = 3Lme(rs) + e(Lmga)- 
If we multiply (3.11) by R,./3 and subtract from (3.9), we obtain 


(3.12) Ait) @ dy?) + de + = W - = olen es 
Applying the properties of the norm and using conditions (i) and (ii), we get 
fe SAY TH+ Welra) | + PW eCLm Pr) Il + 2 [em a) | 


Rn 


+ Mie FI ebm Poll +S" IJ M+ I em au) I 


<5 (AL + Bn) 10M + A om + MB I e(Lm gu) | Le 
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Continuing in this way, we have 
1 TE ines 
moll Sze + 10M + he IL, 


lm || S$R.10°M, 
| W || S$ 6(10°)M + R, 10°M, 


hL 272 3 4 
pods (1 eee MEY iy) + Re dda aw) | 
+ $10°M + || wll + ll ell, 
W eluSsa) Il SW eCys*) + LA fo +2 |W + Welle) I, 


272 373 474 
Vetus) Ws eat) Wt (at AE 4 BE 4 BEY ey | 


», wf. Ser. ; 
+ hR, (7 + 6S + i) | (Le qj) | , 
a, 4 3 eke & 7h 
< a = < é ~— 
+ (SR. + $n) 10°M + (F + 3 ie +X) 10°M 





h’ a“) " h’R ) 2 
— + ——~} 10 °M =) 10 °M. 
+ (4 + 12 0 + r 0 

Now, to estimate the rounding error incurred in one integration step, say from 
j toj + 1, we assume that all quantities obtained at the jth step are exact. Thus, 
in (3.13) we set e(y;*) = 0, e(Lmgs;) = 0, and e(y;) = 0. Denoting the one-step 
rounding error by «6(yj+:), we have 


4 3 di ie rt 252, 
| a(yju) || S ( Rn + 3 h) 10°M + (5 + 54 hRn + a) 10°M 
(3.14) A 
h’Rn 
4 


hee ae 
+ (K+ Zee.) 10 M+ 


= 


10° °M. 


Since h = 10°", we can take R,, = h and get 


13 1 


73 ge + io |, 


PS el 


’ 1 41 —s—m —s—2m 
(3:15) atyha) ls [$10 + 10 + 
which proves the theorem. 
Coro.iary. A bound for the accumulated rounding error, under the hypotheses 
of Theorem 2, is given by 
h) ~~. 
(3.16) Il Cys") Il S ley") le + A - gy Sr) 10°M, 


4 
where f(h) = o (130 + 160h + 100k’ + 27h° + )- 
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Proof. From the definition of y;* we obtain 


ll Cys) Il S ll eCys*) I + =I (Lm gas) || - 
From (3.11), we have 


=e 


Rm || (Lm qj) || S = (10°)M + — 10M, 


< 3, 
2 
Using (3.13), we get 


(3.17) ll e(yFss) |] S el] €(ys*) || + (ih) + fo( Rm, h))10°M, 


where 


fi(h) = 5 (66 + 110h + 64h? + h’), 


folh) = fe (64+ 8R, + 42h + 36h* + 26h' + - ): 


Setting R,, = h, and solving the difference equation corresponding to (3.17), 
we obtain (3.16). 

Remark 6. Theorem 2 gives an upper bound on the one-step rounding error. 
A somewhat better result can be obtained from a statistical estimate of this error, 
if one is willing to make certain assumptions. If we assume (1) that components 
of the errors, e(r;) and e(ZP;) are independent and have a uniform distribution 
between —10°/2 and 10 °/2, and (2) that the bias which would be introduced in 
e(Lqi) and e(Lq) by the coefficient, 4, is eliminated by ‘rounding up’ Lq and 
‘rounding down’ Lq@ , then a direct computation with (3.12) yields as the approxi- 
mate standard deviation of a component of e(y,;*), 


2 5 
(3.18) o; = A E (R + J > (afv/ous)) | 7 -. 
6.3 4% 


This is approximately the standard deviation of an error which is uniformly dis- 
tributed between +R,,10°/2, so that the accuracy is the same as would be ob- 
tained with s + m digits. 

Remark 7. As an example, we follow Gill [1] and integrate y’ = y from t = 0 
to t = 1, with h = 0.1 and s = 6. The results are given in Table 1. The values in 
parentheses are those obtained by Gill’s method [1]. For ¢ = 1, after ten steps, 
we should have the value of e/10. If we use y — g/3 for this value, we obtain 
0.27182810, which is in error by —8 X 10°. (Note that in computing q, the result 
of multiplying by } was rounded up, while in the computation of q , it was rounded 
down. ) 

Remark 8. It is of interest to compare the accumulated error of the above ex- 
ample with the bounds given by (3.16) of the corollary and by statistical esti- 
mates. Since «(yo*) = 0 in the example, and ¢ = Aj = 1, (3.16) becomes 


A f(h) 7 Sth) X 10° 
le(y;*) | S (e Vary x's ae 5) a 
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Comparison of Gill’s Method with the Modified Runge-Kutta Me 


TABLE 1 




















t | Stage r | 2 
| 
| 
5 000 100 000 
he, 250 105 000 
Practing 5 275 105 250 
prslvag 110 525 
| if (110 517) 
ee 8 { 110 517 
a 5 526 116 043 
a. 276 116 319 
| 3 5 830 122 149 
| (122 140) 
02 |; 4 9 { 122 140 
= 6 108 128 248 
i ae 304 | 128 552 
| 4g 6 443 | 134 995 
aera | {134 986) 
0.3 4 9 134 986 
1 6 749 141 735 
2 338 142 073 
3 7 121 149 194 
(149 182) 
0.4 4 12 +t 
7 461 156 643 
2 371 157 014 
3 7369 | 164 883 
: (164 872) 
0.5 4 i a a 
1 8 244 173 116 
2 412 173 528 
3 8 697 182 225 
(182 212) 
0.6 4 13 a one 
1 9 110 191 322 
2 456 191 778 
| 3 9612 | 201 390 
— . | {(201 375) 
0.7 | 4 15 | 01 375 
ar 10070 | 211 445 
ee 502 | 211 947 
3 10 623 | 222 570 
. | { 222 554 
ia ‘ wun 292 554 
bsitanifl 11 128 | © 233 682 
Scat 556 | 234 238 
pee 11 740 245 978 
(245 960) 
0.9 4 18 $m 
1 12 299 258 259 
2 614 258 873 
3 12 974 271 847 
(271 828) 
1.0 4 19 oir eee 
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64 
263 
134 
141 

71 
291 
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78 
322 


164 
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517 


043 
297 
744 


140 


248 
428 
851 
986 
735 
206 
606 
182 
643 
693 
269 
872 
116 
970 
165 
212 
322 
117 
624 
375 
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Now, 





(h) = f0h)_ 130 + 160h + 100h' + 27h! 
Wn) = Teh — 1) 481 + h/2 + W/6 + W/ 24)’ 
and g(.1) = 2.92. Hence, 


| e(y;*) | < 4.98 x 10°. 


To obtain a statistical estimate, we might assume that accumulated error is 
the sum of the one-step errors and that these errors are independent. Using (3.18), 
the standard deviation for one step is about 3.4 X 10°. The standard deviation 
after ten steps is 1/10 times this, or about 1.1 X 107’. 

It is of interest to compare the above results with those obtained from the classi- 
cal Runge-Kutta method, (1.5)—(1.9). These are tabulated below. 


t y 





.100 000 
110 517 
122 140 
134 986 
149 183 
164 873 
182 213 

201 377 

| 222 556 

| 





245 963 
271 831 


moooooocoocoso 
SiON OU RwWNHHO 


4. Floating-point Arithmetic. Since many modern automatic computers provide 
“floating-point” operations, and since the finite-precision modification must be 
applied in a slightly different way when floating-point numbers are used, it seems 
worthwhile to devote a short section to this subject. 

Let us begin by establishing certain conventions. A “digital number in normal 
floating-point form” consists of two parts, a modulus and an exponent. The modulus 
is an aggregate of s decimal digits, the decimal point being placed at the extreme 
left and the digit in position one being non-zero. The exponent consists of two digits 
and represents the power of ten which multiplies the modulus. An algebraic sign 
is associated with each modulus and exponent. Thus, the fixed-point number, 
.00113, would be written as +.113-02 in normal floating form, and —11.3 would 
be written as —.113 + 02. In floating-point arithmetic some of the shift operations 
of formulas (3.1)—(3.5) will be carried out automatically by the positioning which 
must take place in the process of addition or subtraction. The rounding error will 
then be governed by the magnitude of A and the relative magnitudes of y and y’. 
If hy;’ < y;, then Theorem 2 will apply to the procedure (4.1)-(4.5) below, it 
being understood that all errors must be considered as relative errors. If hy,’ > y; 
for some j, the theorem no longer holds. Nevertheless, over an interval, there should 
be a preponderance of points for which hy,’ < y;, so that (4.1)-(4.5) should 
reduce the overall rounding error. 

To explain the meaning of the symbols L and R in (4.1)—(4.5), we must first 
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point out that in automatic computers, the exponent of a floating-point number is 
placed to the left of the modulus. Thus, a shift right m places will not affect the ex- 
ponent if m < s. Now, in the computation of r; , the exponent, x, of the quantity 
in square brackets is compared with the exponent, p, of z;,. If u < p, then R® = 
R,., and L® = L,_,. If u = p, then R® = Ry and L® = Ly. With this inter- 
pretation of the shift operations, the finite-precision modification for floating- 


point arithmetic is as follows: 


20 = Yi; 
(4.1) qo = Qi, 
Po = hf (20), 
r= LR BP» ” ql, 
4 = 2% + ‘1, 
(4.2) a= 3rn — ($Po in qo), 
P, = hf(a), 
= L°R®B(P: — w)], 
= 2 + le, 


(4.4) 


| 
(r2 
(4.3) 2 
P, 
| 


(4.5) | Yin 


94,541 


—t2 — 3 + $Pi, 
hf(z2) — PQ 
LR [P,, 

“a+ 1s, 
—tz+@, 

hf(zs) + 2P2, 


LRP, + al, 
a=2a+%, 


3irs — ($Ps + @s)]. 
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A Statistical Study of Randomness Among the 
First 10,000 Digits of 7 


By R. K. Pathria 


1. Introduction. In connection with the application of Monte Carlo methods to 
various problems in mathematical physics and the drawing of random samples in 
statistics there arose a demand for the so-called random digits. As a result of the 
rapid progress made in these fields of investigation this demand has increased 
considerably during recent years. Consequently, a number of standard sets of such 
digits have been produced and are being put to frequent use by workers engaged in 
these fields {1|-(4]. 

At this stage it appears worthwhile to investigate, as has recently been sug- 
gested by the author [5], the extent to which one can utilize the digits appearing in 
the decimal development of the various constants of mathematical analysis, such 
as e, , etc., for the purposes mentioned above. It is obvious that such a suggestion 
would have been hardly of any practical interest if it had been made at a time when 
the values of these constants were not yet available to a reasonably large number 
of decimal places. However, certain computations of this type have been carried 
out during recent years and in the near future they are to be extended to the point 
where they will surely provide sets of digits as large as the existing ones.* 

Obviously, the question of randomness of the digits to be studied here cannot 
be decided on a priori grounds. One has to subject them to various tests and ob- 
tain internal evidence for their randomness before they can be declared fit for 
practical use. It appears worthwhile to mention here that apart from the specific 
purposes indicated above, a study of this type is fascinating also because of its 
intrinsic interest. It was apparently for this latter reason that Reitwiesner [6], 
at the suggestion of von Neumann, computed the values of z and e to more than 
2,000 decimal places and Metropolis, Reitwiesner and von Neumann [7] carried 
out a statistical treatment thereof by studying the frequency distribution of the 
various digits. This study was extended to about 3,000 decimal places by Gruen- 
berger [8] in the case of e and by Nicholson and Jeenel [9] in the case of 7. 

In the present paper a report is given of the results obtained by applying the 
four classical tests of Kendall and Smith (the frequency test, the serial test, the 
poker test, and the gap test) [10] and a fifth one due to Yule (the five-digit sum 


Received February 27, 1961. 

* Dr. D. B. Gillies of the Digital Computer Laboratory, University of Illinois, has kindly 
informed the author that in a year or so they will probably compute one million digits of e. 
At present their computation extends to 60,000 decimal places. A statistical study of these 
digits is also being carried out by the author and will be reported shortly. 

Very recently Dr. Shanks and Dr. Wrench have made an IBM 7090 calculation of both 
x and e to 100,000 decimal places. The frequency distribution of the decimal digits of both the 
constants has also been computed. The author is highly grateful to Dr. Wrench for illuminat- 
ing communications on this subject. 
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test) [11] to the first 10,000 digits of x.f The value investigated here is the one 
computed by Genuys [13] using the formula 


_<s (-1) { 16 4 \ 
— -A ee 
bao (2k + 1) (S*+* 39-23 H1) ” 

and computing its terms on the IBM 704. The present analysis has been carried 
out mostly in blocks of 1,000 digits each, with a view to discover ‘patches,’ if any, 
that suffer from lack of local randomness. Of course, blocks which are found 
patchy are not suitable for drawing a random sample when used by themselves. 
They have to be suitably diluted by combining them with some of the neighboring 
blocks in order to obtain larger ones which could safely be employed in a statistical 
investigation. 

In comparing the actual frequencies with expectations the x’ test has mostly 
been employed; the rejection levels, following Kendall and Smith [2], have been 
kept at 1 and 99 per cent. 





2. The Frequency Tests. The 10,000 digits of x — 3 have been divided into ten 
consecutive blocks of 1,000 digits each and the frequencies f; with which the various 
digits i(= 0, 1, --- , 9) appear in these blocks have been recorded. These frequen- 
cies, along with the respective values of the statistic x’ and the corresponding prob- 
abilities P for nine degrees of freedom, are given in Table 1. It is only in the case 
of the third and the ninth blocks that the value of P is found to be significant; in 
the former case the deviations from the expected frequencies are too high, while in 
the latter they are too low. 

Taking the table as a whole, of the 100 frequencies recorded 34 deviate from 
the expected value of 100 by more than the standard deviation «(= +/90 = 9.487) 
and 6 by more than 2c. These figures compare well with the corresponding ones, 
namely, 31.73 and 4.55 per cent, for a normal distribution. Further, in the case of 
total frequencies the x’ value (9.318 for 9 d.f.) may be partitioned into three com- 
ponents, with the following obviously satisfactory results: 

















Classification x? | df. P 
Odd versus even digits 0.360 1 ~55 % 
Within groups of odd digits 4.358 4 | -~85% 
Within groups of even digits 4.602 4 ~35% 





3. The Serial Tests. These tests are employed with a view of looking for any 
evidence of serial association among the digits under study. The relevant test here 
consists in classifying the digit pairs (7j) with respect to the members i and j com- 
prising a pair and comparing the frequencies thus obtained with expectations. We 
have tabulated the frequencies for the 10,000 overlapping pairs, formed by the first 





+ Gruenberger [12] has shown how the tests given by Kendall and Smith can be applied to 
any set of digits, punched on IBM ecards, mechanically and without regard to the order of the 
digits on the cards, using standard IBM equipment. In the absence of such a facility, however, 
the author has made the various tabulations by hand and has satisfied himself about their 
correctness by applying suitable cross-checks. 
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TABLE 1 
Frequency Distribution Among the Firsi 10,000 Digits of x — 8 

















\ pisit | 
as a. «6 6 9-8 | # | Pe%) 
Block\ 
a | 
1 | 9 116 103 102 938 97 9 95 101 106| 4.74| ~85 
2 | 8 96 104 86 102 108 106 102 101 106) 4.94| ~85 
3 | 7 gf 9% 7 1 MO 102 9% 108 120/ 22-80] <1 
4 |103 120 10 108 8 102 % 9 95 99| 7.58| ~60 
5 |104 103 88 91 103 108 115 111 87 | 9.38| ~40 
6 | 9 4 98 113 10 97 10 118 90 88| 9.28| ~40 
7 |100 107 98 114 89 108 89 88 98 109) 7.84| ~55 
8 | 97 100 119 95 107 104 108 92 84  94| 8.80| ~45 
9 |101 103 100 103 101 99 98 97 90 108| 1.98) >99 
10 |113 90 110 90 102 113 107 87 94 94) 9.32| S40 
(1-10)* | 968 1026 1021 974 1012 1046 1021 970 948 1014| 9.318) ~40 





* The cumulative frequencies obtaining in this row are in complete agreement with the 
ones given by Dr. Wrench (private communication). See also J. W. Wrench, Jr., ‘The evolu- 


tion of extended decimal approximations to z,’’ The Math. Teacher, v. 53, 1960, p. 644-650; 
v. 55, 1962, p. 129-130. 


10,001 digits of x, in Table 2. The following relations exist among these frequencies: 
2» fis = N 
and 


2d Sim = DX Snn + €m 5 


where N = 10,000 and e¢,, which represents the “end effects’’ is equal to zero if the 

digit m appears either at both the ends of the set or at none; it is equal to —1 if the 

set opens with m and +1 if the set closes with m. In the case under study, we have 

€é; = —1 and « = +1. As a final check on the entries in this table, one verifies 

that the sum Zz fis(t — 7), which should obviously be equal to the difference be- 
1.3 


tween the first and the last digits of the set, is really equal to —5. 
Now, the overall expectation m,; of f:; is, for each of the pairs, equal to Np’, 


where p is the probability of occurrence of a particular digit. The variance of f;; 
is, however, given by 


= Np*(1 + 2ps6.; — 3p’), 
where 6;; is the Kronecker delta. Thus, whereas the expectation for each of the 
hundred elements of the array is 100 on the basis of perfect randomness, the stand- 
ard deviation for the diagonal elements is 10.82 and that for the non-diagonal ones 
is 9.85. The observed values of the root-mean-square deviation are 9.76 and 8.78, 
respectively. Comparing the differences with the standard error in the dispersion 
one finds that none of these values is significant. 

Several essentially equivalent values of x° have been computed from Table 2. 
First, assuming all the hundred types of pairs to be equally likely (expected value 
of 100 for each cell), a x’ of 78.84 is obtained which, for 90 d-f., is at about 80 per 
cent probability level. Second, given the row sums and assuming the ten digits to 
be equally likely to follow (e.g., expected value of 96.8 for each of the cells in the 
first row), a x’ of 69.39 is obtained which, for 90 d.f., is at about 95 per cent prob- 
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TaBLe 2 
Frequency Distribution Among the First 10,000 Overlapping Pairs (ij) of x(= 3.14 --- 78) 











MA 0 1 2 3 4 5 6 7 . 9 | Total 
| | 
0 | 85 103 98 103 98 89 101 93 83 115 | 968 
1 | 99 99 103 102 121 95 106 90 98 113 | 1026 
2 | 101 115 110 99 82 118 100 46101 = 100 95 | 1021 
3 | 102 92 86 94 114 100 90 ©6102 97 98 | 975 
4 | 9% 100 100 89 102 110 103 «(«108~Ss«101 104 | 1012 
5 | 9 117 110 96 108 96 1145 107 96 109 | 1046 
6 | 107 95 117 97 101 124 91 101 90 98 | 1021 
7 | 9 105 99 91 92 101 95 97 103 98 | 970 
8 | 86 97 99 93 96 106 114 83 80 93 | 947 
9 | 112 103 9 110 98 107 106 88 100 91 | 1014 





Total | 968 1026 1021 974 1012 1046 1021 970 948 1014 | 10000 





ability level. Third, assuming the expectation of a particular cell to be one-tenth 
of the corresponding column sum, we get x’ = 69.26 which, again for 90 d-f., gives 
P =~ 95 per cent. Fourth, fitting all the expectations to both the row sum and the 
column sum, a value of 59.83 results which, for 81 d-f., is at about 96 per cent prob- 
ability level. All these figures are obviously satisfactory. 

Next, we have computed from Table 2 the value of the quantity ij whose theo- 
retical expectation and standard deviation are given by 


E(ij) = ()° 
and 
o(ij) = {v — (t)*}-N™”. 

The actual value of this quantity turns out to be 20.062 which deviates from the 
expectation by an amount —1.1 times the 8.D. The probability of equal or greater 
divergence of either sign is about 27 per cent—a result that is not significant. 

So far we have been discussing the question of serial association between the 
neighboring digits comprising the whole set of 10,000 digits. We shall now study 
the various blocks one by one and see if they are individually also locally random. 


For this purpose, we give below the results of the x’ test, carried out on the as- 
sumption of equal a priori probability for each of the hundred cells: 





son oe Bawt ne eas ce Te ae ee 10 





+e 96.6, 82.2) 115.2) 101.4] 96.2) 135.4) 90.8 93.0 80.6} 100.2 


P(%) 30 | 71 + 20 | 3il 0.1; 46 40 | 75 22 


The P value in the case of the sixth block is too low and leads to its rejection out- 
right.{ The only other block for which the P value is rather low is the third one; 
this, however, has already failed to meet the frequency test. 


t It may be noted that this block passed the frequency test very well. The failure here is 
mainly due to an essentially non-random arrangement of the digits in the block. For instance, 
the pair (77) appears 28 times (including 2 triplets and 3 quartets). Such an extreme pattern 
is dangerous even if diluted by one of its neighboring blocks. It can only be made harmless 
by combining with many other blocks. 
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4. The Poker Test. The 10,000 digits of + are printed in 2,000 hands of five 
digits each. Among these hands we have noted the frequency of occurrence of those 
hands whose digits, with respect to their values, are either (i) all different, or (ii) 
one pair and the other three different, or (iii) two pairs and one different, or (iv) 
one triplet and two different, or (v) one triplet and one pair, or (vi) one quartet 
and one different, or finally (vii) one quintet. As usual, the frequencies thus ob- 
tained are compared with expectations. The results are given in Table 3. None of 
the x? values is found to be significant. 

An interesting observation may, however, be made here. Since the deviations 
in the third and the fourth blocks are, on the whole, in the same direction, a group- 
ing of these two consecutive blocks results in a P value of about 1.5 per cent which 
is pretty low, though not below the rejection level. If, however, the rejection level 
were at 5 per cent, as might be the case in a more serious applicatior of these 
digits, this combined sample of 2,000 digits would no longer be considered locally 
random. In that case it would be essential to combine this sample with one of a 
sufficiently large strength before one could employ its digits in an investigation. 


5. The Gap Test. Next, a frequency count has been made of the lengths of the 
gaps between successive zeros of the set. This frequency distribution is compared 
with the expected one only for the whole set and not for the individual blocks, be- 
cause the frequencies in the latter case are too small unless, of course, a very coarse 
grouping of the classes is adopted. 


TABLE 4 
Length Distribution Among the Gaps Between Successive Zeros of x 

















' 


| 
Length of | Actual Expected | Length of the Actual Expected 
the Gap | Frequency | Frequency | Gap Frequency Frequency 
o | 85 | 100.00 16 18 18.53 
1 | 78 | 90.00 17 13 16.68 
> | = | a 4 18 14 15.01 
3 68 | 72.90 | 19 13 13.51 
4 | go | 65.61 | 20 10 12.16 
5 | 61 | 59.05 | 21 13 10.94 
6 47. | 53.14 22 10 9.85 
7 49 | 47.83 23 9 8.86 
8 38 | 43.05 24 7 7.98 
9 42 | 38.74 25-29 32! | 29.39 
10 | 2 ‘| 34.87 30-34 ¢ | 17.36 
11 | | 6-380 | 31.38 35-39 | 8 10.24 
Ss i | fa 4 40-49 13¢ 9.64 
a 25.42 | = 50 | 115 5.15 
1 =. | an 2 - —— 
15 14 20.59 | Total 967 1000.00 





16, 10, 3, 7 and 6, in order of increasing length. 

23, 3, 0, 1 and 2, in order of increasing length. 

31, 2, 3, 2 and 0, in order of increasing length. 

42 each, of lengths 40, 41, 43, 44, 46 and 47; 1 of length 48. 

5 Actual lengths are 51, 51, 54, 55, 60, 62, 63, 65, 65, 65 and 67. 








194 R. K. PATHRIA 


We have 968 zeros in our set and hence 967 gaps; their length distribution§, 
along with the one expected on the basis of perfect randomness, is given in Table 4. 
The x’ value for the grouping as indicated in the table is 28.06 which, for 30 degrees 
of freedom, gives P ~ 55 per cent—a result that is not significant. 

The mean length of a gap (excluding those of length zerc) is found to be 10.190 
which deviates from the expected value of 10 by 0.601 times the standard deviation 
(viz., 0.316). The result is obviously satisfactory. 


6. The Five-Digit Sum Test. This test, as applied here, consists in taking the 
sum of the five digits comprising a (poker) hand as the variable, denoted by the 
symbol i, say, and comparing its distribution over the various hands with the one 
expected theoretically. The latter may be obtained in an elegant manner as follows 
(refer to the alternative approach of Yule [11)). 

If 2; denotes the number of ways in which the five digits of the hand can give a 
sum 7, it will be enumerated by the generating function 


¥ ox = [D2] 


i=0 k=0 


'S. 2)5.(1 —2)* 


lI 


f(x), say. 
It immediately follows that 


- - (5\ (i — 10r + 4 
a= Bc )C“* 4), 


> 2; = Ltf(z) = 10°, 
7 zl 


Moreover, 


being the total number of ways in which a hand of five digits can be formed out of 
the digits of ten kinds. The probability p; for the value 7 of the variable is then given 
by 

Q; 


ae 


which leads to the expected distribution. This is the same as the one given in Table 
I of reference [11]. The mean value of 7 is 22.5 and its standard deviation is 





10 °Q; , 


(41.25)"” = 6.4226. 


The standard error of the mean of n observations is, therefore, equal to 6.4226 
n ”. Further, the standard error of the standard deviation turns out to be 


(18.1/n)"* = 4.2544 n~”, 


§ As a check, it has been verified that the total length of the 967 gaps, as tabulated here, 
is 8988 which, together with the 31 digits preceding the first zero and the 13 digits following 
the last one, makes 9,032—the number of non-zeros in the set. 
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TABLE 5 
Five-Digit Sum Distribution Among the First 2,000 Hands of x 
hye 
Actual Frequencies in 

Expected Fre- Actual Fre- | Expected Fre- 
i quency in a Block Blocks of 400 Hands Each’ quency in the quency in the 

of 400 Hands | I | I m1 | IV | Vv | Whole Set | Whole Set 

Bak Hlth ies ks aa 
0 0.004 | oH 08 8 O AO 0 0.02 
1 0.020 000 ao oO 0 0.10 
2 0.060 o 1 OF O O 1 0.30 
3 0.14400 =| OF OF OF O O 0 0.70 
4 0.280 | Oo oo 1 oO Be 4 1.40 
5 0.504 | 1) O O 1| O th 2.52 
6 0.840 | OF 2 2 oO I = 4.20 
7 Loe. .1.. 2.83 etd 9 6.60 
8 1.980 1 5 te 3 11 | 9.90 
9 2.860 5 Oo} 6] 6 4 ll 14.30 
10 3.984 2 3 i wt Qi 9 | 19.92 
11 5.360 7| dl 9 8 4! 35 | 26.80 
12 6.980 3 7210 6 OB 34.90 
13 8.820 1) 7 9 9 13 49 | 44.10 
14 10.840 | 8 13) 12) 12) 7 ss | 54.20 
15 12.984 | 21) 12| 17) 17) 17 $C 64.92 
16 15.180 13} 10} 16, 18) 15) ees 75.90 
17 17.340 13} 15] 19) 17| 12 76 CO 86.70 
18 19.360 17} 17| 25| 23] 23 105 96.80 
19 21.120 21| 24 18| 26) 24! 113 Cs 105.60 
20 | 22.524 17} 23; 20) 19) 19) 98 | 112.62 
21 | 23.500 21; 19} 23) 22) 32 117 | ~~ 117.50 
3 | 24.000 | 29) 29) 28 20) 29) 135 | 120.00 
23 | 24.000 23| 22| 25| 281 23} 121 | 120.00 
24 | 23.500 | 20) 17| 24| 26) 21) 108 =| = 117.50 
25 22.524 24} 29) 24) 20) 22) 119 | 112.62 
26 | 21.120 | 26, 19) 19} 15, 19) 98 | 105.60 
7 | 19.360 | 19 26 15) 21) 22 103 96.80 
28 | 17.340 | 18 15) 15) 23] 13) 84 86.7 
29 15.180 | 20 14] 10) 12) 12 68 75.90 
30 12.984 | 12 13] 15) 11) 9 60. 64.92 
31 10.840 | 10, 17} 10) 11) 12) 60 54.20 
32 8.820 13} 10) 17) 4 7| od 44.10 
33 | 6.99 | 5 7 4 7 9 32 | 34.90 
34 5.360 667] (6 CU} 8 32 26 .80 
35 3.984 | 5 5} 3 8} 2 18 | 19.92 
36 2.80 | 4 1) 3 3 | 12 | 14.30 
37 1990 | 2 2 4 3 2 13 9.90 
38 1.320 | o oO i ai 3. | 6.60 
39 0.840 | oOo} 2 Of 0 1 ?_ 4.20 
40 | 0.504 | 11 Oo OF OF O a 2.52 
a4 0.280 0 0 Oo oO 1 1 1.40 
42 | 0.1440 | OF 08 OF OF 0 0 0.70 
43 | 0.060 | 0 O OF OF O 0 0.30 
44 | 0.0200 | 0 OF 0 O O 0 0.10 
45 | 0.004 | 0 OF O 0 O 0 0.02 
2000 


Total 400 .000 400) 400) 400, 400 400 
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TABLE 6 
Mean Values of the Sum i, with Differences from Expectation, Etc. 


| 














| j 
. Difference from | Divided by | Square of the 
Block | Mean Value | Expectation | Standard Error 'Preceding Column 
! | | - 

I | 22.7450 | +0.2450 | +0.7629 | 0.5820 

II | 22.7550 | +0.2550 | +0.7941 | 0.6306 
Ill | 22.4625 | —0.0375 | -0.1168 | 0.0136 
IV | 22.0925 —0.4075 | -—1.2690 | 1.6104 
Vv | 22.1800 —0.3200 | —0.9965 0.9930 

The whole set | 22.4470 | —0.0530 | —0.3690 
TABLE 7 


Standard Deviations of the Sum i, with Differences from Expectation, Etc. 

















| | 
Block | Standard (Difference from| Divided by | Square of the 
| Deviation | Expectation | Standard Error 'Preceding Column 
ae |— | | 
I | 6.3945 | —0.0281 | —0.1321 | 0.0174 
II 6.4475 | +0.0249 | +0.1169 | 0.0137 
Ill | 6.2919 | -0.1307 | —0.6143 0.3773 
IV | 6.2241 | -—0.1986 | —0.9334 0.8713 
V | 6.3716 | —0.0510 | —0.2398 | 0.0575 
The whole set | 6.3524 | -0.0702 | -0.7379 | 





The actual frequency distribution obtained from the 2,000 hands of the set is 
given in Table 5, where the results are also given for consecutive blocks of 400 
hands each, i.e., comprising 2,000 digits each. The actual distribution is compared 
with the expected one through the mean values of the variable and its dispersions. 
In Table 6 we have listed for each of the five blocks, I to V, and for the whole set, 
the mean values and their deviations from the expectation in terms of the standard 
errors of the mean. None of the various deviations is found to be significant. In 
fact, the chance of equal or greater divergence, of either sign, in the case of the 
whole set is about 70 per cent. Moreover, even if we group together the last three 
blocks (each having a deviation of the same sign) the corresponding result comes 
out to be about 17 per cent. Still worse, if we take the last two blocks, for which 
the deviations are not only of the same sign but also of the greatest magnitude, the 
result is still about 11 per cent. Further, we note that the sum of the entries in the 
last column of the table is 3.83. Entering the x’ table with this value of x’ and 5 
degrees of freedom, we find P to be about 60 per cent. 

Finally, we study the standard deviations in the value of the variable as ob- 
tained from the frequencies tabulated above and compare them with the corre- 
sponding theoretical expectations. The relevant figures are given in Table 7. Ex- 
pressing the deviations in terms of the standard errors of the standard deviation, 
we obtain results which do not exceed unity. Further, entering the x° table with 
the sum of the squares of these numbers, namely, 1.34, and 5 degrees of freedom, 
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we find that P lies between 90 and 95 per cent. For the whole set, the deviation 
of the actual standard deviation from the expected value is equal to —0.74 times 
the corresponding standard error; the chance of an equal or greater deviation of 
either sign is about 46 per cent. 
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University of Delhi 
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An Extended Table of Roots of 
In'(%) Yn'(8%) — Jn'(B%) Yn'(x) = 0 


By John F. Bridge and Stanley W. Angrist 


An eigen-equation that frequently occurs in mathematical physics involving 
an annular cavity is 


(1) Jn'(%)¥n' (Ba) — Jn'(Bx)¥,'(x) = 0 


where J,(z) and Y,(z) are respectively Bessel functions of the first and second 
kinds. 

J. McMahon [1] gave an asymptotic expression for the roots to this equation, 
D. O. North [2] obtained a root smaller than the first root given by the asymptotic 
expression of McMahon, and R. Truell [2] developed a graphical method for ob- 
taining this root. 

H. B. Dwight [3] gave the first six roots of equation (1) for values of n from 1 to 
3 and for various values of 8 from 1 to 4. 

The purpose of this paper is to extend the range and accuracy of the roots to 
equation (1) and to determine the ranges of the solutions for which the asymp- 
totic expression proposed by McMahon is sufficiently accurate. 

The calculation of the roots was accomplished by trial and error substitution 
in the following equation: 


(2) F(x) = Jn’ (x) Yn’ (Bx) — Jn’ (Bx) Y,'(z). 


Starting with x = 0.1, F(x) was calculated, increasing x in steps of 0.1 until 
F,(2) changed sign. A linear interpolation was then used to determine the ap- 
proximate value of the root, then the Newton-Raphson iteration procedure was 
used until two successive approximations of the root value were within +10°°. 
The root thus obtained was compared with the root obtained with all four terms of 
J. McMahon’s asymptotic expression (3) 





2 3 

: So, 2-8 . * — e+ 
af ao4 & + EE 4 2 ++: 

o te 2) _m+3 

B-1’ a 
_ 4(m* + 46m — 63)(6° — 1) 

3(88)*(8 — 1) 
." 32(m® + 185m” — 2053m + 1899)(6° — 1) 


5(88)(8 — 1) recs 
If the two values were within +2 X 10°°, the asymptotic value was used for 
that root and all larger roots for the given value of n. If not, the procedure was 
continued until the next root was found and compared with its asymptotic value. 


6 m = 4n? 





(3) 
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EXTENDED TABLE OF ROOTS OF Ja'(z)Y¥2'(62) — J,'(6x)Y¥,'(x) = 0 
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The roots calculated by the McMahon expression (3) are indicated in the tables 
by the letter A printed after the number. It should be noted that equation (3) 
does not give a root for s = 1. 

The great speed of the IBM 704 digital computer was used to advantage in the 
solution of this problem. The Bessel functions were generated by using recursion 
relationships. Starting with an arbitrarily small number for a Bessel function of 
high order for a given argument, successively smaller orders were calculated until 
Bessel functions of 6- to 7-place accuracy were obtained. It was not possible with 
the particular procedure used to evaluate functions with arguments larger than 50. 

The authors wish to express their appreciation for the advice given by Pro- 
fessor L. 8. Han who originally suggested the problem. 
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Polynomial and Continued-Fraction Approxi- 
mations for Logarithmic Functions 


By Kurt Spielberg 


1. Introduction. In this article we present the coefficients of approximations 
which are well suited for the calculation of logarithms on digital computers. The 
approximations have been derived by means of the IBM 704 program IB CTR. 
They are chosen so as to approximately minimize the absolute error over the appro- 
priate interval of the argument. The method is described in detail in references 
[1], [2]. 

Similar selected polynomial approximations have been made available by C. 
Hastings [3]. The approximations of the present article, however, cover a much 
wider range of accuracy and should allow the coding of efficient double-precision 
subroutines. 

Continued fraction approximations have been used systematically by E. G. 
Kogbetliantz and the author in connection with subroutines for the IBM 704 and 
709 computers (see e.g. [4], [5], [6], which contain many references to other litera- 
ture on rational approximations). The reader should note that the continued 
fraction approximations given in this paper not only allow for computation with 
fewer second-order arithmetic operations (multiplications and divisions) but also 
are intrinsically more accurate than polynomial approximations with equal num- 
bers of constants. 


2. Polynomial Approximations. In the case of digital computers, the argument 
can be assumed to be in normalized floating point form: 
A. Binary machines: 


(1.a) y = 2*f 
i --- integer, f --+ fraction, (3) S$f<1. 


B. Decimal machines: 


(1.b) y = 10'-F 
I --- integer, F --- fraction, (js) = F <1. 
The natural logarithm is then evaluated in accordance with the relations: 
(2.a) log. y = (t + log: f) -log. 2 
(2.b) log. y = (1 + logy F) -log, 10. 


To obtain efficient polynomial approximations, one starts with the well known 
series 





(3) log. © + 


U 


* = 2A(x/v) + (2°/30*) + (2°/50°) + +++] 


- | 
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which converges in the interval [—v < x < »]. Since we intend to “economize” 
the power series by means of Chebyshev polynomials, we focus our attention on 
the interval [—1 < x < 1). To this end we introduce the rational transformations 


Se ee) 


von (Vite a {Ea 
(4b) F =U (7+), z=aV-( 


and determine u,v, U, V so that the interval [—1 < zx < 1] maps one-to-one onto 


the intervals [(4) < f S 1] and [(+s) S F <S 1] respectively. 
The parameters are determined from the endpoint conditions: 


v= (V2 +1) 


_ (V+1\_,, (V-1\_1 — 
(5b) U (Ft) -1,0 (Fa4)" 2 wat Rie ns 
V = (V/10 + 1)*/9. 


On substituting these values into equation 3, we obtain the following power series 


for logef and logiF : 
(6.a) logs f = 2-log, e-[(x/v) + (2°/3v°) + ---] — (4) 
(6.b) logi F = 2-logye-[(x/V) + (2°/3V*) + ---] — (3), [(-lsza lJ 


The 704 program IB CTR is now applied to produce polynomial approximations 


to the functions log, f + (4) and logy F + (4). These approximations have the 
form: 











m 
-(m) _2i—1 
> Cei-1°L 


i=l 


_ ,f— (vV2/2) 1 
(7) x ” =r (9/2) (73/2) for ioe: S + (5) 


_ pF —(v/10/10) ( 
x YF + (10/10) for logi F + ;). 


For computational purposes, however, it is preferable to introduce the variables 


zx x 
-Orz=—s,; 
v i 


Sn* (x) 





z= 


(8) WP hidd ul D oft. 2 


In Tables 1 and 2 we give the coefficients c§”’ for those m which result in approxi- 
mations of less than or equal to 16-digit accuracy. IB CTR performs operations 
to 16-digit accuracy only. Its primary output, however, consists of the increments 
Aa; which have to be added to the power series coefficients a; of the given function 
to produce the coefficients of the approximation polynomial. We therefore give 
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tables of these increments from which the reader can construct approximations of 
great accuracy by simple hand computation. 

All approximations of the form (7) have been tested at more than 100 points 
in the interval [—1 < x < 1]. Instead of the complete error curves we submit, 
for simplicity, three “error parameters.” 

E, --- a theoretical upper bound of the magnitude of the absolute error caused 
by a truncation of a Chebyshev series to m terms 


E, --- the maximum magnitude of the absolute error encountered in the de- 
scribed test 
0 
E;--- >> aoi-1, the maximum absolute error incurred by a truncation of the 
i=m+1 


given power series to m terms. 


The sets of increments have been tested as follows. From the definitions we 
infer that (for z = 1) 


> Qin + p> Qin = >» (@ei1 + Aae:1) + max (E,, Ep) 


i=m+1 i=l 


or 


E; = yi Aaa; + max (E,, Ez). 

Selected tests of this type have consistently been satisfactory. The reader should 
note, however, that these tests do not usually apply to the last two digits due to 
the unfortunate fact that EZ; has been printed only to 6 digits. In order to obtain a 
better check, at least up to “triple precision accuracy” on the IBM 704 (2~”), 
we have therefore coded a triple precision logarithm subroutine based on the given 
increments. The accuracy of the subroutine was verified by an application to 
functional relationships of the form log (x-y) = log x + log y. We have every 
reason to believe that all of the given increments will be found to be completely 
accurate. 


3. Continued Fraction Approximations. An approximation polynomial can be 
transformed into a rational approximation with the same number of constants by 
means of the “multiple truncation procedure” described in [2] and implemented 
in IB CTR. It is shown in [2] that the rational approximation may actually be con- 
siderably better than the original polynomial approximation. The results sub- 
mitted in the present article furnish an excellent instance of this behavior. 

Rational approximations can readily be transformed into continued fractions 
which can be evaluated in fewer operations. In Tables 3 and 4 we give the con- 
tinued fraction expressions for (log, f + 4) and (logis F + 4) up to 16-digit accur- 
acy. They are of the form 





. * Gi | ae & © = Gimj2_| 
Om (2)/e= Ho + Boa + araet + FF Hea’ 
where m = 3, 4, --- and [3m] is the largest integer < = . For even m, the con- 


stant Ho is zero. 
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All continued fractions have been checked at more than 100 points in the inter- 
val [—1 S z = 1]. These checks were executed in double precision arithmetic, 
i.e., with wordlengths of 16 digits. The user of a particular continued fraction 
approximation should briefly analyze how much round-off error may accrue on his 
machine due to limited wordlength and subtraction of numbers of equal magnitude. 
For large H; and G; serious loss of accuracy might occur in this manner. In the case 
of the logarithmic functions, however, little difficulty should arise from this source. 


4. Use of Tables. To illustrate the use of the tables we give a few simple ex- 
amples. 
a) Polynomial approximation, three coefficients: 


/2 
—cieceei ue ja i 
~Sjs z 
> oe =— *) 


loge f = fa*(z) — } + (.32) 10 = exe + ge + ese” — 4 + (.32) 10” 
Table 1: C; = 2.88539 12843 --- 
96147 14921 --- 
Cs = .59895 53187 --- 


Another way of writing the approximation would be 


C3 


1 z-—l1 
— < 2 Zzg= 
5 Sts V2, 5) 


logex = cyz + c92° + cg” + (.32) 10” 





b) Continued fraction approximation, three coefficients: 





/2 

1 io 

aD 4 f & 1, = 

2 /2 

f+> 

loge f = gs (z) — Ste (.48) 10° = z| Hy +=s4 Gh an ae 48) .10° 
2 2+ A 7 aa 

Table 3: Hy = 1.29200 70987 
H, = —1.65676 26301 
G, = —2.63985 77031. 


c) Use of increments: 
The increments of Tables 1’ and 2’ should be added to the coefficients: 
(2 logs e, $ logs e, 2 loge e, ----) 
and 
(2 logio e, $ logy e, logis e, ----) 
respectively. 
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In our computations we have used the constants 
2 logs e = 2.88539 00817 77926 8146 
2 logwe = .86858 89638 06503 6553. 
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TECHNICAL NOTES AND SHORT PAPERS 


Some Remarks on Modular Arithmetic and 
Parallel Computation 


By H. S. Shapiro 


1. Introduction. A question that has been discussed in recent years is that of 
parallel computation. Can a given computation be broken up into independent 
assignments which may be performed simultaneously? Traditional methods of 
computation are almost entirely serial, with the consequence that one cannot con- 
vert extra computing capacity into significantly greater speed. Thus far the only 
general method which has been proposed for achieving parailelism is the use of 
“modular arithmetic’”—that is, for some collection of relatively prime integers 
m, , -** ™ one performs the calculations (mod m;) independently; the final result 
is then obtained by solving a system of simultaneous congruences. Such a procedure 
is possible provided that (i) the calculation consists entirely of additions and mul- 
tiplications of integers* (sc that the corresponding calculations (mod m;) are jus- 
tified) , and (ii) each number sought in the calculation is an integer known a priori 
to lie in an interval of length < mm, --- m,. 

Modular arithmetic, when applicable, has the advantage of being free from round- 
off errors; moreover, addition and multiplication (mod m) are carry-free. Another 
feature is that in some types of calculation (for instance, tabulation of the values 
of a polynomial for equally spaced values of the argument) the calculation (mod m) 
is much simplified by the periodic repetition of the values being calculated. It there- 
fore seems of interest to show how computations of practical importance may be 
carried out within the limitations (i) and (ii) above. In this note we discuss division, 
linear equations, and the first boundary value problem from the standpoint of 
modular arithmetic. 


2. Division. Let us consider the problem of finding d binary digits of the quo- 
tient 3 where x and y are integers, 0 < x < y. The most natural approach is to 


choose an n > 0, and let r be the least non-negative residue of 2" (mod y), then, 





writing N = , adh , N is an integer easily computed (mod m) if (m,y) = 1, and 
we have 

tN <% 

ay’ 


these numbers differing by % < x . Hence, the integer xN, converted to binary 


Received May 22, 1961. This report represents results obtained at the Institute of Mathe- 
matical Sciences, New York University, under the sponsorship of the Office of Naval Research, 
United States Navy. 


* Division is allowable only when the modulus is relatively prime to the denominator, and 
the quotient is an integer. 
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notation, gives d digits of the quotient = , providing 2”** > z. Here the approxi- 
mation is from below; by working instead with 
N’ = 2+8 
y 
we get an approximation from above. The objection to this procedure is that cal- 
culation of the residues of N(mod m;,) is simple only in case (m;,y) = 1, and so 
each denominator gives rise to a certain set of “forbidden” moduli. 


The following division algorithm is free from this defect. Let b = 2" be the 
smallest power of 2 not less than y. Then the sequence ¢, defined by 


binnr = (D-— ye +z 


’ SMF, 


converges to : (for an arbitrary choice of t). Writing s, = b"t, , we get 


(1) Saat = (b — y)8, + "2. 


If so is chosen to be an integer, all the s, are integers, easily computed and periodic 
(mod m), for all m without exception. In place of (1), we can use the convergent 
iteration 


(2) Sati = —(y — @)8,’ + az 


where a is the greatest power of 2 not exceeding y. By choosing the better of (1), 
o~9a1-—¢ 
b a 
numbers is <4) we achieve good convergence. The necessary a priori estimate of 
8, (or 8,’) respectively, and the degree of approximation after n iterations, are 


readily obtained, and from this the magnitude of M = Im; sufficient for calcula- 


(2), i.e., (1) or (2) according as 








is smaller (and at least one of these 


tion of ~ to the required accuracy is known. Since ¢, arises from s, upon division 

by 2™, i.e., shifting of a binary point, the conversion of s, from modular to binary 
x 

notation gives the initial digits of = directly. 


A variation of this division algorithm which gives a simpler recurrence at the 
expense of an auxiliary calculation is this. Suppose that the (given) binary expan- 
sions of z, y are 


xz = a2” + a2” + --- a, 

y = bo" + 2°" + --- db, 
where a; , b; are 0 or 1. We may suppose a, = b, = 1 since multiplication and divi- 
sion by powers of 2 is trivial. Then » = 2” “f(4), where f denotes the function 
Gy + dpat + --- aot? 
by + bert + --- dott 


It is easy to show that, in the Taylor expansion (3), we have | c, | < nth Fibonacci 


(3) f(t) = =l+qtitetf+-:--. 
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number, so that (3) converges at least* for | t| < vd Frat = .618 ---. Moreover, 


from (3), the integers c, are readily computed in terms of the a; , b; by a recurrence 
obtained upon cross-multiplying in (3). 
But in terms of the c, the division can be carried out, using the scheme 


(4) Sai = 2 Sn — Ca - 


If so is an integer (say, s) = 0) the s, are integers and 


teil 1 ea 
new 2 (5) y 
once again, the a priori bounds on s, , and the rate of convergence (which is uni- 


formly rapid with respect to all divisions) is readily obtained. 
These algorithms may be adapted, in an obvious way, to any radix. 


3. Linear Equations. Given the system of k equations (in vector notation) 
(1) Az = b, A = || a; || 


where the a;; and 6; are assumed to be integers, the direct adaptation for modular 
arithmetic is to compute d = det A, and replace (1) by the system 


(2) Ay = db 


for the integer variables y; . Operating with moduli m; such that (d, m;) = 1 (as- 
suming, of course, d ~ 0) the solution of (2) (mod m,) is very simple (say, by 
Gaussian elimination). Crude a priori bounds on the y; may be obtained (e.g., by 
Hadamard’s determinant inequality) when they are not available from physical 
or other considerations. However, there is again the objection that this scheme 
allows “forbidden moduli” which vary with the given problem. Moreover, the 
solutions x; are found as quotients, necessitating divisions which are non-trivial. 
These difficulties disappear when an iterative method is employed. Let us suppose, 
to keep the discussion simple, that by preliminary transformations (1) has been 
put into a form where | ai | > > >;2:| ai; |. Let ec; be the integer defined by: (i) 
c; has the same sign as aj; ; (ii) | c; | is the least power of 2 not less than | a;; |. Let 
C be the diagonal matrix of the c; , and c = max | c; |. Then the system (1) may 
be rewritten in the form 


(3) Cz = Dx + b 
(writing D = C — A), and the iterative scheme 
(4) Crna. = Da, +b 


converges to the solution of (1), because of the supposed diagonal-dominance of 
A. Finally, putting y, = c"x, we get 


(5) Yngr = (€C*)Dyn + (cC)be". 


* The fact that the series (3) converges at least in this circle, i.e., that the polynomial in 
the denominator cannot vanish in this circle, was discovered in another connection and pointed 
out to the author by D. J. Newman. 
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Since cC has integer entries, all components of all y,, are integers (if ye is so chosen) . 
The iterative scheme (5) is suitable for modular computation; we omit a detailed 
discussion of rate of convergence and a priori bounds. 


4. First Boundary-Value Problem. When the equation V*u = 0, subject to given 
boundary conditions, is solved numerically by the method of (square) nets, it is 
customary to use an iterative method of solution which is known to converge at a 
rate that can be estimated in terms of the geometry of the region. This problem is 
in principle subsumed in the above discussion, but two factors make it especially 
simple from the standpoint of modular computation. First, the transformation to 
integer variables is particularly simple and especially favorable to a solution in 
binary notation (owing to the special significance of the number 4 in this iteration). 
Second, the maximum principle gives a good a priori bound on the solutions. For, 
writing the iterative scheme symbolically as 


Ungar = Aun + D 
and setting 4"u, = vn, we get 
(1) Ung. = 4Av, + 4"*b. 
Suppose that the components of b are integers (i.e., that the given boundary values 
are integers, which is achieved by shifting a binary point). If, then, the initial values 


Vo = Uo are chosen to be integers lying between the smallest and greatest boundary 
values, all components of the v, computed from (1) are integers, lying in the range 


4"B, =< v < 4"B, 


where B, and B, are the min (and max) of the prescribed boundary values. 


5. Remarks on Other Iterative Methods. There are many other important itera- 
tive methods in numerical analysis, but not all of these seem well adapted to modu- 
lar computation, because in many cases transformation to integer variables leads 
to integers that are too large to be computed practically, i.e., an excessively great 
number of moduli are required. We may illustrate this with a simple example. 
Suppose we try to solve the equation 


(1) e+zr—%=0 
by means of the convergent iteration 
In41 = —2n +3, % = 0. 
Letting 2", = y, we have 
(2) Yan = 2 x, yo = 0. 


The numbers y, defined by (2) are integers, whose initial binary digits coincide 
with those of the positive root of (1). Moreover, calculation (mod m) of the num- 
bers y, from (2) is quite trivial. However, since y, is of the order of 2” even 10 
iterations using moduli of the order of 50 would involve us (roughly) in calculating 
a 1000 binary digit number, by solving a system of 200 simultaneous congruences— 
a lot of work to solve a quadratic equation. Newton’s method would lead to the 
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same difficulties. A feasible method for solving quadratic equations by modular 
computation can, however, be based upon the Taylor expansion (1 — 42)” = 
Zana )e”. 


6. Concluding Remarks. Preliminary analysis indicates that parallel computa- 
tion, using modular arithmetic, is feasible for certain kinds of problems. The parallel 
computation envisioned here leads very swiftly to a solution encoded “in modular 
notation.”’ By this is meant, a system of simultaneous congruences, whose solution 
(in a specified interval), written as a binary number, has as its initial digits the 
binary number which is the goal of the computation. For results of practical value 
it will probably be necessary, at the very least, to use moduli whose product exceeds 
10°°. Hence the feasibility of rapid solution of large-scale systems of congruences 
will determine the timesaving possibilities of the method. Any a priori knowledge 
about the solution, such as might be obtainable from a preliminary rough solution, 
analog computation, etc., leads to a reduction in the number of necessary moduli, 
i.e., knowledge of r binary places reduces the product of the m; needed by a factor 
2”. Again, in such a case as the boundary value problem, where the values of the 
solution at neighboring net points differ by amounts which can be bounded a priori, 
this fact might lead to a considerable reduction of labor in the ‘“‘conversion” phase 
of the problem. 
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Permutations with Restricted Position 


By Frank Harary 


In his book on combinatorial analysis, Riordan [4, p. 163-164] discusses permu- 
tations with restricted position and mentions an open question: 
“Any restrictions of position may be represented on a square, with the elements 
to be permuted as column heads and the positions as row heads, by putting a 
cross at a row-column intersection to mark a restriction. For example, for per- 
mutations of four (distinct) elements, the arrays of restrictions for the rencontres 
and reduced ménage problems mentioned above are 
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1 2 3 4 1 2 3 4 
1 x 1 x x 
2 x 2 x x 
3 x 3 x x 
4 x 4 x x 
recontres ménages 


Since a square of side n has n° cells, and a cross may or may not appear in each 
cell, it is clear that with n elements a problems are possible (this includes per- 
mutations without restriction, for which no cell has a cross). However, many of 
these are not distinct since, from the enumeration standpoint, the relative rather 
than the absolute position of the crosses is important; for example, all n’ prob- 
lems having just one cross on the board are alike. The exact number of distinct 
problems, for any n, is not known, but some progress in this direction will appear 
in this chapter.” 

In this note, we show that the question has been virtually solved in [2], and 
shall obtain an explicit formula for P, , “the exact number of distinct problems, 
for any n.” For we shall see that the chromatically nonisomorphic bicolored graphs 
with n points of each color, which are enumerated in [2], are in a one-to-one cor- 
respondence with the distinct problems involving permutations on n objects with 
restricted position 

A binary matriz is one in which every entry is 0 or 1. Consider the set M of all 
square binary matrices of order n. We say that two matrices A and B in M are 
equivalent if B can be obtained from A by the following three operations: 

1. Perform any permutation on the rows of A, obtaining A, 

2. Perform any permutation on the columns of A; , obtaining A» 

3. Either leave the matrix A, as it stands or take its transpose, obtaining 

A; = B. 

Obviously, this is an equivalence relation and it is clear that these three opera- 
tions are independent. This equivalence relation partitions M into equivalence 
classes. The number of distinct problems of permutations on n objects with restricted 
position is thus the number of equivalence classes of the matrices in M. In the above 
quotation from Riordan [4], the presence of an z in his matrix corresponds to that 
of a 1 in the associated binary matrix, while a blank space in his matrix becomes a 
0 in the binary matrix. 

A graph consists of a finite collection of points together with lines joining certain 
pairs of distinct points. When two points are joined by a line, they are adjacent. 
A graph is said to be colored with k colors if each point is assigned one of these 
colors, any two adjacent points have different colors, and all & colors are used. A 
bicolored graph is one which has been colored with two colors. 

To a given restricted permutation problem represented by the binary matrix 
A = (a;;), there corresponds the bicolored graph with 2n points 1, 2, ---, n, 1’, 
2’, ---, n’ in which point 7 is joined by a line to point 7’ if and only if a;; = 1. Thus 
for n = 4, the rencontres and reduced ménage problems give the bicolored graphs 
of Figure 1. 
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4—— 1 4 
2——- 2 2 
?... 9 3 
4-—-~ 4’ 4 

recontres menages 
Fie. 1 


Two graphs are isomorphic if there is a one-to-one correspondence between 
their set of points which preserves adjacency. Two bicolored graphs are chromatically 
isomorphic if there is an isomorphism between them which preserves color, i.e., 
two points of the first graph have different colors if and only if their image points 
do. Clearly, chromatic isomorphism is an equivalence relation on bicolored graphs. 
As an illustration, we show in Figure 2 all bicolored graphs (up to chromatic iso- 
morphism) in which there are two points of each color, together with the corre- 
sponding matrices. 

Lemma. Two square binary matrices are equivalent if and only if the corresponding 
bicolored graphs are chromatically isomorphic. 

Proof. We translate the three defining operations for equivalence of matrices 
into graphical terms. Any permutation of the rows of a matrix A in M corresponds 
to a renumbering of the 7 points of the first color in the associated bicolored graph 
G. A permutation of the columns of A becomes a renumbering of the n points of 
the second color. Finally, transposing the matrix A amounts to interchanging the 
two colors assigned to the points of G. Clearly, these three operations serve to char- 
acterize chromatic isomorphism. 

A formula for the counting polynomial 


(1) Gan(Z) = p> box" 


where b, is the number of chromatically nonisomorphic bicolored graphs with n 
points of each color and g lines which have been found in [2]. Let P,, be the number 
of inequivalent matrices in M, i.e., the number of distinct types of restricted 
permutation problems (on n objects) ; then 


(2) Pa = gnn(1) = bo +b, + es + dae. 


For example, we see from Figure 2 that goo(z) = 1 + x + 2x° + 2° + 2‘; hence 
P, = 6. The number P,, may be found from the cycle index of the ‘“‘exponentiation 
group” S;? (where S, is the symmetric group of degree n) using the enumeration 
lemma of [1]. This is the same procedure as substituting 1 for x in the formula for 
9nn(x), which is derived using Pélya’s method [3]. To give the result conveniently, 
we require the following notation: 
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(t) = (t1,%, ---, i.) denotes a partition of n such that: 


(3) i, + 2g + --- + ni, = 2, 
(4) »(i) = |] ki! 
k=l 


and d(r, s) is the greatest common divisor of the positive integers r and s. 
1. ¥ 1 7 1 1 1 Y 


i 1’ 
= ta) ee eS 2 2 o Sue es * 


1 


ae 
b db db dE dE 











1 
. a 
0 0) 
Fic. 2 


We may now state the formula for P, . Let 


n 


(5) a(i,j) = 2) tjad(r, 8) 

(6) (3) LE; + ()] +> +144 2) | + ¥ iried(r, 8). 
Then 

(7) Pp, =! > 1 gat. 4 1 _ 1 oa 





~ 3 hk W@rGQ) 2 & (3) ~ 


where the first sum is taken over all pairs (7), (j) of partitions of n and the second 
sum is over all partitions (7) of n. 

We illustrate for n = 3 whose three partitions ; , m2 , 7; are 1 + 1 + 1,1 + 2, 
and 3. These may be written as the sequences 


(3,0,0),  (1,1,0), (0,0, 1) 
respectively. The values of a(a;, 7;) for n = 3 are given in the matrix: 


{ 


alm. 7;) 


~ 

wa © 
i a) re ey 
wns w 
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while the values of »(x;) and 8(2;) are given in the table: 





i v(4:) B(x) 
1 6 6 
2 2 2 
3 3 2 


Hence we have 


P=! aie are er Ph ee aod 
=-|— = a —+ — = ~|— = — | = 26. 
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Evaluation of the Zeros of Cross-Product 
Bessel Functions 


By L. Jackson Laslett and William Lewish 


1. Introduction. There is considerable interest in the zeros of certain cross- 
product Bessel functions which arise in solving Bessel’s equation subject to Dirichlet 
or Neumann boundary conditions at r = a, b, 


(la) J,(qa)Y,(qb) — Jn(qb)Y,(qa) = 0 
or 
(1b) Jn'(ga) Yn'(qgb) — Jn'(qb) Y,'(qa) = 0, 


because of their well-known application in physical or engineering problems for 
which the use of cylindrical coordinates is appropriate. In many instances attention 
may be directed primarily to the zeros of such functions when n is not large because 
of the interest in the lower-order modes which are possible in the physical problem 
under consideration, but cases may also arise in which the higher-order modes will 
warrant attention in order to determine the circumstances in which such possibly 
unwanted modes may become excited. 

Solutions to (la) and (1b) have been discussed by a number of writers [1], [6], 
and results presented in the form of algebraic formulas, in tables, or graphically. 
For application to problems in which (b — a)/(b + a) is small and in which n may 
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be large, however, it appeared appropriate to make an independent investigation 
of the initial roots of (la) and (1b) by study of characteristic solutions of Bessel’s 
equation in the interval a < r < 6 without explicit reference to the usual Bessel 
and Neumann functions. Approximate analytic formulas have been obtained from 
which estimates may be made of the characteristic values, for the case of the first 
Dirichlet root and for the first two roots subject to the Neumann boundary con- 
dition, and an independent numerical determination of the characteristic values 
and characteristic functions has been made with the CYCLONE electronic digital 
computer at Iowa State University for cases in which (b — a)/(b + a) was given 
the values 0.001, 0.01, and 0.1. It is the purpose of the present note to summarize 
the results of this investigation, for which more detailed results will be available 
elsewhere (see Section 5). 


2. Transformation of Bessel’s Equation. It may be noted that, due to the nature 
of the customary Bessel functions of high order, and in particular because the func- 
tion J, remains quite small until its argument is comparable to its order, the lowest 
characteristic values, g, will be in the neighborhood of n/b for n large. For this 
reason, and to focus attention on the interval a <= r < b, it is convenient to define 





b-—a 
(2a) ied = 
(2b) 5=n | (« b+ey = “|, 
and 
(2c) o ot — (b+ a)/2 


b-—a 
In terms of these quantities, 


(3) r= oF (1 + a2), with —-lsz 





lA 


1, 
and Bessel’s equation assumes the form 


(4) flat) 2] + [aca + we) +7 + ™ et Z.gint-z |Z = 0. 
dx dx 


1+ a 

The solutions to (4) which are of interest are those for which the Dirichlet 
boundary condition (Z = 0) or, alternatively, the Neumann boundary condition 
(dZ/dx = 0) applies at = +1. When the Dirichlet boundary condition is applied, 
it may be convenient for some purposes to make the transformation 





(5) S = (1 + az)""Z, 
in terms of which (4) may be written 

d's (9°/4) + 'n(2 + wt|s Ay 
(8) dat * [ id (i + a2) * 


with S(+1) = 0. 








228 LASLETT AND LEWISH 


Physically, it is seen that the quantity 7 which is introduced here denotes the 
ratio of the width (b — a) to the mean diameter (b + a) of an annular region. 
For 7 only slightly less than unity, the annular region extends substantially from 


r =0 to r= b and the roots g . as o of (1a) or (1b) may then be expected to be- 


come one-half the corresponding roots, u, of the simpler equations J,(u) = 0 or 
J,'(u) = 0, respectively. 

For 7 < 1, the terms in (4) or (6) which contain 7, save in some cases those 
which involve the combination 7°n’, may either be ignored in determining simple 
analytic formulas for 6 or may be treated as a perturbation. 





3. Approximate Analytic Formulas. For 7 < 1, the characteristic values, 6, for 
(4) or (6) may be obtained by a perturbation method [7] in which the unperturbed 


equation is taken as simple harmonic, provided n is not too large. In this way we 
find 


7a) For the first Neumann root: 6 ~-=7'n’ — 


The nature of the characteristic solution associated with the first Neumann root is 
such that it is very nearly constant when 7°n’ is small. In such cases the form of 
the solution is approximately given by Z ~ 1 + 7°n? (: _ =) Similarly, the first 
Dirichlet and second Neumann solutions are respectively of the general character 
cos 5 x or sin x. The region of applicability of (7a—-c) may be considered to be that 


for which 7°n’ < 1; of equal or greater interest, however, are the results for the 
case n’n’ > 1, which is discussed below. 
In cases for which y’n’ is not small, but 7 < 1, it may suffice to replace (4) by 


2 
(8) ei + [6 + 2n'n’-2]Z = 0. 


Solutions of this approximate equation may be written explicitly in terms of Bessel 
and Neumann functions of order }. It then follows, moreover, that for 7°n° at least 
somewhat greater than unity (e.g., n’n” > 6) the solution of interest is substantially 


(,w2 g on : 
ig | J (<5) + Ji (£5) |. for & 
(9) Z 3n'n' 3n’n 
‘ oe a 3/2 
es it lg Wea (: | & | ) ‘ for ¢ 


3n'n? 


IV 


0, 





IIA 


0, 


where ~ denotes 6 + 2n’n"z, since the first Hankel function then becomes sufficiently 
small at s = —1 as to satisfy adequately the boundary condition normally imposed 
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at that point. The characteristic values, 5, may then be estimated by application 


of the desired boundary condition at z = 1, aided by tables of Js:;; and J2/s {8}, 
(9], 


(10a) For the first Neumann root: 6 ~ —2n'n® + 1.617249'n*”, 
(10b) For the first Dirichlet root: 8 ~~ —2y'n? + 3.711519 n*’, 
(10e) For the second Neumann root: § ~ —2n'n’ + 5.156199'n“*. 


The numerical constants which appear in (10a, b) are seen to be, as expected, twice 
the numerical coefficients given in series developments for the first maximum and 
first zero of J,, when n is large [9 (Sect. 15.83, p. 521)]; [10 (Sect. VIII.3.6, p. 143)]. 
Characteristic values for solutions to (8) must necessarily be somewhat less nega- 
tive than —2n’n’ in order that the coefficient of Z be positive for some values of x 
in the interval —1 < x < 1. For 7'n’ large, the characteristic solutions are rela- 


tively large only for values of x near unity, in a region whose width is roughly two 
or three times (n°n”)~*. 


4. Computational Results. The differential equation (4), suitably scaled, was 
integrated with the CYCLONE digital computer at Iowa State University, using 
the Runge-Kutta process [11], [12]. Runs were made for several values of n, with 
7 given in turn the values 0.001, 0.01, and 0.1. In each case the value of 6 was ad- 
justed, by trial, to give solutions satisfying the desired Dirichlet or Neumann bound- 
ary conditions. A larger number of integration steps was employed to traverse the 
interval —1 < x < 1 in cases in which n° was large, since more rapid changes 
of the function occur in certain portions of that interval in such cases. The effect 
of truncation error was found, by tests in which the interval size was halved, to be 
sufficiently small that use of the finer interval only affected the final value for the 
function or its derivative (in the Dirichlet or Neumann cases, respectively) by less 
than 10° of the maximum value and the consequent error in 6 could thus be judged 
when tabulating the results of the investigation. 

The characteristic values 6 determined computationally are listed in Table I. 
By comparing calculated values of 6 obtained for (7a-c) and (10a-c) with the 
values in Table 1, the accuracy of (7a—c) and (10a-c) can be ascertained. See Table 
VI [13]. Figure 1 depicts the nature of the associated characteristic functions, for 
n = 0.01, for various representative values of n. Since the contribution from 6 makes 
a relatively small change in the characteristic value for the original Bessel equation 
when n is large, use of (2b) in connection with the values of 6 given in Table I should 
afford accurate characteristic values for q in such cases. In the application to phys- 
ical problems it is interesting to note from Figure 1 the features mentioned in Sec- 
tion 3, namely that at small n the first Neumann solution does not show a pro- 
nounced variation with x and the other characteristic solutions have approximately 
the form of circular functions, while at large n the characteristic solutions become 
large only in a small interval near x = 1. 


5. Availability of Detailed Results. The analytic work of Section 3 is presented 
in greater detail, and the computational results reproduced directly from the tele- 
printer output of the CYCLONE, in an Ames Laboratory report [13]. The report 
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1 
also includes approximate values of / Z'dx, suitably normalized with respect. to 
=’ 


the value of Z or dZ/dx at x = 0 and at x = 1, for 7 = 0.0001 and for representa- 
tive values of y°n’ in the range 0 through 20. This report is available from the Office 
of Technical Services, U. 8. Department of Commerce. Two copies of the report 
have been deposited in the file of Unpublished Mathematical Tables which is main- 
tained by Mathematics of Computation and may be made available on loan to in- 
terested individuals. 
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On the Computation of Lommel’s Functions 
of Two Variables 


By J. Boersma 


In 1942 Zernike and Nijboer [1], [2] introduced a new expansion of Lommel’s 
functions of two variables in connection with calculating the diffraction integral 
of a circular aperture. In this article it is shown that this expansion is very well 
suited for the computation of these functions. (The author is much indebted to 
Dr. Bottema of the Physical Laboratory of the University of Groningen, who drew 
his attention to this formula.) 

Lommel’s functions of two variables are defined in the following way (Cf. [3], 
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formulas 16.5 (5) and (6), p. 537, 538), 


Ca) v+2m 
U,(w, z) > (-0"(¥) Jr42m(2) 


(1) 


2 
V,(w, z) = cos ¢ + 5 + 3) + U_.42(w, z). 


The present article deals with the computation of Lommel’s functions of two 
variables of integral order. Owing to the recurrence formulas, (Cf. [3], formulas 
16.5 (7) and (8), p. 538), 


ikea) + Ceaahiglie (“) Iz) 
(2) 


V,(w, z) + Vi42(w, z) - (“) J_.(z), 
it is sufficient to compute Lommel’s functions for two successive integral values of v. 


The first table of Lommel’s functions of two variables of integral order is to be 
found in Lommel’s memoir on diffraction at a circular aperture [4]. Lommel gives 


2 2 2 2 
tables for Uw, 2), se U2(w, z), and for ~ Vo(w, z), A Vi(w, z) to six decimal 


places for values of the arguments w = x(x)10x, z = 0(1)12, and w = x(x)8z, 
z = 0(1)12 respectively. 

Quite recently, a table [5] by Dekanosidze has been published which gives 
tables of Ui(w, z), U2(w, z), Vi(w, z), Vo2(w, z) to six decimal places for a some- 
what uncommon domain of values of the arguments: 


w = 0.5(0.02)1.2(0.05)4(0.1)6.2, z= w(0.01)4V/w 
w = 6.3(0.1)10, z= w(0.01)10. 


The tables may also be used outside this domain of values by means of the rela- 
tions (Cf. [5], formulas (7) and (8) ) 


U.(w, 2) = (—1)"Ve (Z,2) 
w 


Va(w,2) = (—1)"Un (E,:). 
Ww 


Dekanosidze’s tables have been computed by means of a power series expansion 
2 


(3) 


of U,(w, z) and V,(w, z) in powers of 5s (Cf. [5], formula (3)), 


U,(w, z) = > (—i)" — = (< y Uy4m(w, 0) 


m=0 2w 
(4) 2\m 
= zZ duce bs 
V,(w,z) = > (— +6) V,—-m(w, 0). 


When » is integral, the coefficients of these power series contain a factor of 
the type U,(w,0) and V,(w,0), where n is an integer. U,(w,0) and V,(w, 0) are 
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given by [5], formula (4), which contains some printing errors. (Cf. [3], formulas (Cc 
16.52 (11)—(16), p. 540). The correct formulas are as follows: 








w\™ 
_ n < m (5) 
Usa(w,0) = (—1) | e085 - 2X (- 1) () | hei 
- a | 
Uenyi(w,0) = (—1)”" P -= (~ ” Gm +1! 
. w. nr 
5) U_,(w, 0) = cos ( + _ " 
Vo(w,0) = 1, Vass(w,0) = 0 lie “ 
| ve 
a tir 
-_ fu 
V_2n(w, 0) — —1)" > (— (2m )! al 
ee 
=_ | cé 
V_2n-1(w, 0) aj (—1)” > (— Sy + 1)!" { v 


The computation of expressions (5) for not too small values of w may suffer a 
from loss of digits owing to the alternating character of the series. This same ob- e: 
jection arises in computing the alternating series (4) when z is not small. d 

We now turn to Zernike’s method. Here Lommel’s functions of two variables W 
are expanded in products of Bessel functions [ 


1 MA 
Ui(w, z) + iU2(w, z) = wei | Jo(zt)e mt dt 0 
(6) e 


we'“* 5 i"(2n + 1) zs La (5). 


The essential advantage of this expansion is the following. In (6) all terms of the 
infinite sum have an absolute value smaller than 1 for all real values of w and z, 
so, contrary to Dekanosidze’s method, there is no danger of loss of digits. This is 


readily proved by applying the recurrence formula for Bessel functions (Cf. [3], } 
formula 3.2(1), p. 45) L 


2(2n + 1) Senses) = Jon(z) + J on42(z) 
hence 
l 
is 
Da | 


Similarly for J,4;(x) the following integral PORE TI is valid: 


[(2n + 1) 2m @) Jan (2) | + =| Jonge (2) | € 1. 


Ly 
2 | 


> et 
Insy (2) = (—8)" 4/5 J eP,, (t) dt, 





COMPUTATION OF LOMMEL’S FUNCTIONS OF TWO VARIABLES 235 


(Cf. [3], formula 3.32(2), p. 50) which may be estimated by 


\Jun (2)| 52 4/Z = 4/%; 
Vera) V5 eo} 


Another advantage is that U; and U2 are calculated simultaneously because 
each of them is found by adding alternate terms of one single expansion. 

The Bessel functions of odd and semi-odd order which are required in equation 
(6) may be computed very suitably by means of the recurrence technique de- 
veloped by Goldstein and Thaler [6]. When computing a table of Lommel’s func- 
tions on an electronic computer, it is possible to store these sequences of Bessel 
functions, after which various values of w and z may be combined to give U;(w, z) 
and U,(w, z). 

The method may still be used, even for large values of w and z, though in that 
case a rather large sequence of Bessel functions must be computed. 

A comparison of the two methods has been made for the case w = 20, z = 20. 
When Dekanosidze’s method was followed for both functions U;(w,z) and U2(w,z), 
a total of 31 terms of the series in equation (4) had to be taken into account, 
each term being computed to twelve digits in order to obtain an accuracy of four 
decimal places (hence a loss of eight digits). When the method described here 
was followed, 11 terms were already sufficient to give simultaneously U,(w, z) and 

U2(w, z) with the same accuracy without any loss of digits. 

Finally, the method described here has been used to recompute Lommel’s 
original tables [4] (the functions Vo(w, z) and V,(w, z) have been computed by 
equation (1)), the results being given below. In these tables, the decimals which 
deviate from Lommel’s values have been italicized. Besides that, all values of 


hence 


2 V,(w, z) differ by a factor —1 from Lommel’s values because the definition of 


V,.(w, z), as used in the present article and in [3], differs by a factor (—1)” from 
Lommel’s original definition. (See the footnote at the bottom of [3], p. 537.) 
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TABLE OF LOMMEL’S FUNCTIONS OF TWO VARIABLES 





























2 2 Uste, 2) 2 vate, 5) | = Vote, =) | 2 Vi(w, 2) 
0 +0.636620 +0 .636620 | +0.636620 0 
: +0.539802 +0.580638 | +0.479744 | 2. 088772 
+0.298890 +0.433460 | +0.055002 | —0.213022 
3 +0 .032426 +0 .247396 | —0.383136 —0.055402 
5 ~ 0.173625 | ~0.022058 | 0-4s0e40 | Lo-d52883 
—0.173625 | —0.022258 | +0.45 | +0.252583 
6 —0.094496 —0.056474 +0.278235 | —0.636026 
7 +0.011819 —0.041090 | —0.676734 | —0.023425 
9 ee | Tees | touee| feune 
-Uo ° brs é —U. 
10 +0.007050 +0.017334 | +0.148505 | +0.630010 
ll —0.035803 +0.007209 —0.245499 | —0.620117 
12 —0.039518 —0.004302 +0 .504943 | +0.342522 
w= 2x 
| 
P 2 vite, 2) 2 vatw, 2) 2 Velw, 2) | =~ Vilw, 2) 
0 0 +0.636620 | +0.318310 | 0 
1 —0.047572 +0.559947 +0.242644 | —0.022268 
2 —0.156737 +0.362318 | +0.059998 | —0.057118 
3 —0.250135 +0.124194 | -—0.115909 | —0.041158 
4 0.259807 | 0.066375 | —0.159699 | +0.044515 
5 —0.172632 —0.155073 —0.025675 | +0.118189 
6 —0.036806 —0.141277 +0.164916 | +0.050182 
H +0.106459 | 0.007007 | --0-111109 | —0:1s0077 
9 +0.065125 +0.045577 —0.268535 | +0.116651 
S| eee | dome | ieee | Hess 
a 6 . | . e | a 7) 
12 —0.046338 —0.013333 —0. 155667 | —0.331053 
w= 3r 
2 z Ui(w, 2) 2 U2(w, 2) : Vo(w, z) 2 Vilw, 2) 
0 —0.212207 +0.212207 +0.212207 0 
1 —0.221811 +0. 150853 +0.162106 —0.009903 
2 —0.233157 —0.000541 +0.044154 —0.025710 
3 —0.209291 —0.162866 —0.065351 —0.020817 
4 —0.127691 —0.255583 —0.096321 +0.012549 
5 —0.005131 —0.240383 —0.034488 +0.046240 
7 +0.156654 | —o-o1iss | +0.0903%% | 0.025132 
A —0.010133 +0. 45 —0.025132 
8 +0.124279 +0.079105 +0.025841 —0.081135 
9 +0.038751 +0.098281 —0.095896 —0.046848 
10 —0.043777 +0.059355 —0.116648 +0.074775 
11 —0.077056 | +0.001898 +0 .030683 +0.133189 
12 —0.052825 —0.035554 +0.171787 —0.007646 
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TABLE OF LOMMEL’s FUNCTIONS OF Two VARIABLES 
w=A4r 
| | | 
2 = Ui(w, 2) a Us, =) 2 V ow, 2) | 2 Vile, 2) 
0 | 0 | 0 40.159155 | 0 
1 +0.000759 | —0.037360 +0.121669 —0.005572 
2 +0.010698 | —0.122929 +0.034215 —0.014526 
3 +0.043564 | —0.194789 —0.045730 | —0.012219 
4 +0.100101 | —0.196607 —0.068629 | +0.005486 
5 +0.157469 | —0.114657 —0.027959 | +0.024001 
6 +0.179329 | +0.013349 +0.035306 | +0.021696 
7 +0.141278 | +0.122492 | +0.063628 | —0.006592 
8 +0.052246  +0.160931 | +0.029137 | —0.036976 
9 —0.046242 | +0.119653 | —0.038977 | —0.033318 
10 —0.104762 | +0.033666 | —0.072886 | +0.013463 
11 —0.097782 | —0.043578 | —0.027365 | +0.060545 
12 —0.039653 | —0.073719 | +0.061663 +0.044025 
w= dr 
: 2 vs(w, 2) | 2 uste,s) | 2 vetw,s) | 2 Vilw, 2) 
| 
| | 
0 +0.127324 +0.127324 | +0.127324 | 0 
1 | -++0.123693 +0.101421 | +0.097369 | —0.008566 
2 | +0.116978 | +0.043947 | +0.027779 | —0.009316 
3 +0.114163 | +0.000633 | —0.035346 | —0.007972 
4 +0.114193 +0 .008654 —0.053424 +0 .003028 
5 | -+0.103761 +0.068242 —0.022719 +0.014668 
6 | +0.066399 | +0.140571 +0.024568 +0.013915 
7 —0.000892 +0.173623 +0.046307 —0.002302 
8 | —0.077848 +0.137780 +0.024054 —0.020595 
9 —0.128770 +0.046262 —0.021724 —0.021117 
10 —0.125074 | —0.053370 —0.048087 +0.002141 
11 —0.066712 —0.110067 —0.027077 +0 .029850 
12 +0.014393 —0.100916 +0.025355 +0.030738 
w = 6r 
| | ? 
| = Us(w, 2) 2 Us(w, 2) Votw, 2) | = Vile, 2) 
OE, NE a, ee ee ee ee a 
0 | 0 +0.212207 | +0.106103 | 0 
1 —0.005291 | +0.187222 +0.081156 | —0.002477 
2 —0.017713 | +0.128841 +0.023335 | —0.006476 
3 —0.030684 | +0.074204 | —0.028890 —0.005594 
4 | —0.041775 +0.052871 | —0.043819 | +0.001917 
5 | —0.055411 +0.064625 —0.018991 +0.009906 
6 —0.076993 +0.080250 | +0.018958 | +0.009616 
7 —0.103193 | +0.064885 | +0.036479 | —0.000963 
8 —0.118371 +0.006400 | +0.019824 | —0.013120 
9 —0.103083 | —0.072681 —0.014730 —0.014204 
10 —0.050141 —0.129002 —0.035334 | —0.000298 
11 +0.024503 —0.128183 —0.022325 | +0.017291 
12 +0.086795 —0.069153 +0.013475 +0 .020232 
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w= 74 
Z 2 Ui(w, 2) E U2(w, =) z Vow, z) 2 Vilw, z) 
0 —0.090946 +0.090946 | +0.090946 0 
1 —0.092742 | +0.067502 | +0.069570 | —0.001820 
2 | 0.095331 +0.011836 | +0.020097 —0.004761 
3 —0.093184 —0.042950 —0.024469 —0.004136 
4 —0.083669 —0.069547 —0.037187 +0.001326 
5 —0.069496 —0.065233 —0.016277 +0.007149 
6 —0.055115 —0.050886 +0.015516 +0.007029 
7 —0.040559 —0.051440 +0.030185 —0.000452 
8 | —0.019618 —0.073600 +0.016738 —0.009122 
9 +0.014183 | —0.098796 —0.011166 —0.010150 
10 +0.057961 —0.097344 —0.027951 —0.000826 
ll +0.095375 | —0.053091 —0.018474 +0.011275 
12 +0.104159 | +0.020684 | +0.008673 | +0.014010 
w= &8r 
: 2 vxtw, 2) | 2 vate, ») | = Vetw, 2) Vico, 2) 
| 
0 0 0 | +0.079578 0 
1 +0.000190 | —0.018684 +0.060878 —0.001393 
2 +0.002679 —0.061687 | +0.017639 —0.003647 
3 +0.010993 —0.099549 | —0.021243 —0.003179 
4 +0.025878 —0.107904 | —0.032324 +0.000974 
5 +0.043375 —0.084168 —0.014231 +0.005408 
6 +0.057603 —0.046848 +0.013178 +0.005359 
7 +0.065631 | —0.018864 | +0.025804 —0.000228 
S +0.069350 | —0.008873 | +0.014458 —0.006731 
9 +0.071904 | —0.005808 | —0.009042 —0.007608 
10 +0.071829 | +0.009152 | —0.023198 —0.000876 
ll | +0.061297 +0.043412 | —0.015655 +0.007971 
12 | +0.031978 | +0.082866 | +0.006317 +0.010231 
w= Or w= 107 

2 vutw, 2) > vate, 2) 5 tite; ») Us(w, =) 
0 +0.070736 +0.070736 0 +0. 127324 
1 +0.069624 +0.055367 —0.001905 +0. 112361 
2 +0.067676 +0.020712 —0.006385 +0.077695 
3 | +0.067323 —0.007570 —0.011132 +0.046173 
4 | +0.068669 —0.008852 —0.015445 +0.035955 
5 +0.068172 +0.017625 —0.021256 +0.047323 
6 | +0.061099 | +0.053530 —0.031103 +0.063679 
7 | +0.04568/ | +0.076475 —0.044829 +0.065343 
8 | +0.024884 | +0.076750 —0.058321 +0.044755 
9 | +0.003842 | +0.062426 | —0.065889 +0.011066 
10 | —0.014690 | +0.049474 | —0.064357 —0.018761 
11 | —0.032149 | +0.046026 | —0.055052 —0.034099 
12 | —0.050776 +0.044629 —0.041669 —0.037907 




















—_ to Ae 


ADDITIONS TO CUNNINGHAM’S FACTOR TABLE OF n‘* + 1 239 


Additions to Cunningham’s Factor Table of n'+1 


By A. Gloden 


This note is the fulfillment of a plan to present in a readily accessible and con- 
cise form a complete list of additions to the factor tables of n* + 1 published by 
Cunningham [1], which give the prime factors (with certain omissions herein 
supplied) of all such integers not exceeding 1001‘ + 1. Cunningham’s factoriza- 
tions were found with the aid of his tables [1] of solutions of the congruence 


x + 1 = 0 (mod p) 


for p < 10°. 

The subsequent tables of S. Hoppenot [2], A. Delfeld [3], and the writer [4] 
have provided an extension of these congruence tables to include all admissible 
primes between 10° and 10°. 

The factorizations presented in the present note have been extracted from a 
number of sources. The data corresponding to even values of n < 442 and to odd 
values of n < 523 have been published previously by M. Kraitchik [5] and N. G. 
W. H. Beeger [6]. The remaining data have appeared in a series of papers by the 
writer [7]. 

In Cunningham’s table of factors of n* + 1 for n = 2(2)1000 there appear 97 
incomplete entries. Of these, 66 are now identified as primes, corresponding to the 
following values of n: 


320 442 526 616 742 800 952 
328 466 540 624 748 810 962 
334 472 550 628 758 856 966 
340 476 554 656 760 874 986 
352 488 556 690 768 894 992 
364 492 566 702 772 912 996 
374 494 568 710 778 914 
414 498 582 730 786 928 
430 504 584 732 7 930 
436 516 600 738 798 936 


Of the remaining 31 incomplete entries, 14 correspond to primes of the form 
(n* + 1)/17, 


namely, when n = 648, 678, 682, 706, 746, 784, 790, 818, 842, 876, 882, 892, 954, 
988. 

Furthermore, (n* + 1)/41 is a prime when n = 888, 946, and 998. Thus, there 
remain 14 omissions to be considered in Cunningham’s table, for even values of n. 
These factorizations are now given in extenso. 
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426 
59: 

640 
698 
714 
820 
828 
844 
850 
880 
924 
938 
980 
982 


A. GLODEN 


né+1 
129553 - 254209 
203569 - 628193 
174289 - 962609 
189017 - 1255801 
216841 - 1198537 
626929 -721169 
176041 - 2669977 
246289 - 2060273 
170873 - 3054937 
290737 - 2062673 
158993 - 4584689 
809273 - 956569 
780049 - 1182449 
137593 - 6758489 


In the companion table of factors of n* + 1, for n = 
82 incomplete entries, of which 68 have now been shown to correspond to primes 
of the form (n* + 1)/2. The related values of n are herewith listed: 


403 
405 
415 
419 
431 
445 
449 
453 
455 
463 


471 
477 
479 
487 
503 
505 
513 
517 
523 
537 


539 
543 
551 
561 
567 
573 
579 
605 
607 
613 


623 719 
639 721 
643 725 
649 745 
657 761 
677 769 
681 795 
701 805 
703 811 
713 819 


821 
829 
833 
843 
845 
855 
857 
879 
883 
891 


1(2)1001, there appear 


895 
913 
917 
919 
931 
963 
965 
997 


Moreover, (n* + 1)/2-17 is primeforn = 801,859, 865, 869, and 961; (n* + 1)/ 
2-41 is prime for n = 957 and 981. In addition to these entries, it is now known 
that (n* + 1)/2-17° is prime when n = 1001. 

Consequently, there remain only six entries to be considered, and for these the 


complete factorizations of (n* + 1)/2 are as follows: 


n 
565 
595 
685 
889 
893 
941 


(ms + 1)/2 
157217 - 324089 
137321 - 456353 
147377 - 746969 
505777 -617473 


17 - 104233 - 179441 


132961 - 2948521 


In conclusion, I should like to state that this paper was prepared as the result 
of a suggestion made to me by Dr. J. W. Wrench, Jr. that I consolidate my results 


an mm Awe is S65 
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and those of other researchers which complement the factorizations of n* + 1 
published by Cunningham. 


11 Rue Jean Jaurés 
Luxembourg 
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On the Generation of All Possible Stepwise 
Combinations 


By Gary Lotto 


Conventionally, when all possible combinations of all possible subset sizes from 
a set of n are desired, a binary count is performed. Associating the units digit with 
the number 1, the two’s digit with the number 2, the four’s digit with the number 
3, etc., the binary count 0001, 0010, 0011, 0100, 0101, 0110, 0111, 1000, etc., be- 
comes associated with the combinations 1, 2, 12, 3, 13, 23, 123, 4, ete. This is useful 
in such procedures as the analysis of variance. 

The above order of combinations requires that, when computing on data from 
one combination to the next, either (a) the calculation starts anew, or (b) if algo- 
rithms exist for generating a new function from the old one by single steps of either 
including or deleting a number from the combination, more than one step may be 
required. For example, we may go from the combination “2” to the combination 
12” by “including 1.”’ But going from “12” to “3” requires “deleting 1, deleting 2 
and including 3.” 

Given, then, that a problem may be solved for some combination of k elements 
from the solution for the superset of (k + 1) elements or the subset of (k — 1) 
elements, is there an algorithm for generating all possible combinations which goes 
through the fewest recursions? 
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TABLE 1 
i | A(i) B(i) C(i) 
1 | 00001 +1 1 
2 00010 +2 12 
3 00011 os 2 
4 00100 +3 23 
5 00101 +1 123 
6 00110 ~2 1 3 
7 00111 ac 3 
8 01000 +4 34 
9 01001 +1 | 1 34 
10 01010 42 1234 
11 01011 ay 234 
12 | 01100 | ~3 2 4 
13 | 01101 | +1 12 4 
14 | 01110 | ~2 1 4 
15 | 01111 | ay 4 
16 | 10000 | +5 45 





The author has used the following algorithm to generate all combinations of in- 
dependent variables in a multiple regression problem: 

(1) For each step, carry the cycle number i of the combination which is to be 
generated. 

(2) Divide z by 2, then the quotient by 2, etc., until the remainder is not 0. The 
number of divisions performed is k, the number to be included or deleted. 

(3) Divide the quotient of the last division in (2) by 2. If the remainder is 0 
include. If the remainder is 1, delete. 

The algorithm is equivalent to inspecting the lowest non-zero bit in the binary 
representation of 7. If this is the kth bit (counting from the right), the number k 
is to be included or deleted. The (k + 1)st bit instructs inclusion or deletion: if 0, 
include; if 1, delete. 

Define A(z) as the binary representation of 7, B(i) as +k if the number k is to 
be included, or —k if k is to be deleted on cycle 7, and C(7) as the resultant combina- 
tion. Table 1 gives the first 16 values of 7 and these functions. 

Given combinations 1 through (2° — 1), all combinations of (k — 1) elements, 
the additional combinations which must be generated in order to produce all com- 
binations of k elements are reproductions of the first (2° — 1) combinations, to 
each of which has been added the kth element, plus the combination of element k 
alone (in effect, a reproduction of the zero combination, plus element /). 

The algorithm produces these combinations by: (1) including k on the 2*'st 
cycle, and not deleting it before the 2‘th cycle, and (2) reproducing the B(i)’s in 
reverse order with opposite sign (B(2"* + c) = —B(2** — c)), thus on each 
cycle deleting from the combination that which we, 2c cycles before, included into 


b 
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it, or including that which we, 2c cycles before, deleted from it, until the (2°)st 
combination, which corresponds to the empty set plus element k. 

Proof of (1). Since the binary representation of 2°” is a 1 bit followed by (k — 1) 
zeros, the kth element is included on cycle 2°. The kth element will remain until 
the binary number 11 followed by (k — 1) zeros appears. This will be on cycle 
number (2* + 2°") > (2* — 1). Thus, all combinations from 2*~* through (2° —1) 
will include the kth element. 

Proof of (2). Since (2** +c) + (2*"* —c) = 2", the binary representations of 
(2** + c) and (2° — ec) correspond in all their low-order zeros, and the low-order 
1, in which they also correspond. The bit above the 1 must differ in the two num- 
bers, due to the binary carry. Thus, B(2** + c) = —B(2** — e). 

To complete the proof by induction, we may note, by Table 1, that the algorithm 
has generated all combinations for k < 4. 


University of Pittsburgh, and 
American Institute for Research 
Pittsburgh, Pa. 


Generation of Permutations by Addition 


By John R. Howell 


1. Introduction. Suppose one wishes to generate the k! permutations of k dis- 
tinct marks. Representing these k marks by 0, 1, 2, ---, ( — 1) written side by 
side to form the “digits’’ of a base k integer, then the repeated addition of 1 will 
generate integers whose “digits” represent permutations of k marks. Many num- 
bers are also generated which are not permutations. D. H. Lehmer [2] states that 
this so-called addition method can be made more efficient by adding more than 1 
to each successive integer. 


2. Method. In this note, we show that the correct number greater than 1 to 
add to this integer is a multiple of (k — 1) radix k. 

Lemma 1. The arithmetic difference radix k between an integer composed of mutu- 
ally unlike digits and another integer composed of a permutation of the same digits 
is a multiple of (k — 1). 

Considering the process of “casting out nines,” it is obvious that the two in- 
tegers are congruent mod (k — 1). Hence, their difference is zero mod (k — 1). 

The method seems to have two advantages. First, one can generate all k! per- 
mutations in lexicographic order. Second, all permutations “between” two given 
permutations can be obtained. The process can be made to be cyclic if upon ob- 
taining (k — 1), ---, 0 one takes the next permutation to be 0, 1, ---, (kK — 1). 


3. Example. Suppose we wish to generate the 4! permutations of 4 marks. 
Representing these 4 marks by 0, 1, 2 and 3, we add 3 radix 4 to 0123 to get 0132. 
Continuing this process we get the 4! permutations desired. The array below shows 
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the first 16 numbers generated by this process. An asterisk marks each integer 
whose digits represent a required permutation. The other integers were rejected 
because of the occurrence of repeated digits. 


Sequence Integer Sequence Integer 
1 0123* 9 0303 
2 0132* 10 0312* 
3 0201 11 0321* 
4 0210 12 0330 
5 0213* 13 0333 
6 0222 14 1002 
7 0231* 15 1011 
8 0300 16 1020 


4. Adaptation to a Computer. In a computer such as the IBM 7090 where con- 
vert instructions are available it is easy to do radix k arithmetic. Otherwise one 
could simulate the process by adding 9 digit-wise and testing the resulting sum for 
having unique digits each one of which is one of the original k digits. 


5. Acknowledgments. This method was developed when the author was with 
the Statistics Department, Agricultural Experiment Station, University of Florida, 
Gainesville, Florida, in connection with the problem of obtaining a particular ar- 
rangement of the rows of a Latin square. He wishes to thank Mark Robinson of 
Martin Marietta Corp. for suggestions concerning the writing of the manuscript. 


Computer Applications Department 
Martin Marietta Corporation 
Orlando, Florida 
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Multiple Quadrature with Central Differences 
on One Line 


By Herbert E. Salzer 


Abstract. The cocfiicients Adm in the n-fold quadrature formulas for the stepwise 


integration of (1) y = f(z, Y,Y', °°; ww); at intervals of h, namely, for n 
even, (2) 5"yo = h” ms, (1 + Aind”)fo + ---, for n odd, (3) w6"yo = h” Dory 
(1 + Adnd”)fo + +++, are tabulated sanutly for n= 1(1)6, m = 1(1)10. They 


were calculated from the well-known symbolic formulas (4) * ie = (3/D)"f, (5) 


(8/D)" = (dh/2 sinh(8/2))" and (6) w = (1 + 87/4)" =1 +> - sat _ 
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~<8 

a + ---. For calculating y“”, replace n by n — r in (2) and (3). Use of (2) 
and (3) avoids the solution of (1) by simultaneous lower-order systems for n > 1, 
as well as mid-interval tabular arguments, requires only even-order differences, on 
a single line, and provides great accuracy due to rapid decrease of A},, as m in- 
creases. However, the integration may be slowed down by the need to estimate 
and refine iteratively the later values of y, y’, ---, y“” required in 8°"f, . Refer- 
ence to earlier collected formulas of Legendre, Oppolzer, Thiele, Lindow, Salzer, 
Milne and Buckingham, reveals that Thiele and Buckingham come closest to (2), 
(3), as their works contain schemes that involve just tabular arguments through- 
out. For n odd, they give formulas that are based upon the series in &” for 
(1/z)(6/D)” instead of u(é/D)” as in the present arrangement. 


1. Purpose and Scope of Tabulated Formulas. Given a differential equation 


(1) y” _ f(x, y, y’; P y””), 


and a sufficient number of starting values at intervals of h, there are very convenient 
numerical integration formulas for obtaining either 6”yo , for n even, or ud"yo , for 
n odd, in terms of just the even-order central differences of f = f(z, y, y’, ---, y"”) 
at x = 2%, denoted by 8 "fy. This article tabulates the exact values of A}, , the 
coefficients of 6°"fy , for n = 1(1)6, m = 1(1)10, in the following numerical inte- 
gration formulas: 


10 
(2) 5"yo = h” >, (1 + Alnd’”)fo + ---, for n even, and 
m=1 
10 P 
(3) ys"yo = h” >> (1 + Alnd’”)fo +--+, for n odd. 
m=1 


The computation of Az, was based upon the symbolic form of (1), or D"y = f, 
from which 


(4) é"y = (6/D)"f. 
The well-known operational formula, 
(5) (6/D)" = (éh/2 sinh™ (8/2))", 


was used to obtain the coefficients of 6°” in the series for (6/D)". For even n, this 
yielded (2). For odd n, (5) produces integration formulas that express mid-interval 
values of y in terms of tabular values of f. To obtain (3), which involves tabular 
values of both y and f, we multiply (5) by u, giving y, on the left side, the numerical 
interpretation of a mean central operator 4(B'? + FE), and considering yu, on 
the right side, a symbolic even function of 6 according to 


9 


(6) p= (1487/4) =1 += 


8 


a 
128 1024 32768 
Integration of (1) also requires formulas for the stepwise determination of the 

derivatives y“, r = 1(1)n—1. By noting that D”“y” = f, we can still employ 

(2) and (3), as well as the same quantities &°"f, , merely replacing n by n — r. 
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In the use of (2) and (3) we avoid the widespread practice of breaking up a 
higher-order equation into a simultaneous first-order system where each equation 
requires its own set of differences. Also there is no occurrence of formulas involving 
mid-interval arguments. Among the attractive features of this scheme is the em- 
ployment of just alternate even-order differences that are on a single line. Besides 
the concise and economical appearance of (2), (3), the rapid rate of decrease of 
Am With increasing m is seen to provide high accuracy. 

On the dampening side, the user is reminded that the higher-order central differ- 
ences of f(z, y, y’, ---, y“” ) in (2) and (3) involve later values of y, y’, ---, y°” 
that must be estimated at first, probably by some kind of extrapolation. Then (2) 
and (3) might be used in some iterative refining scheme, the details depending 
upon the particular functional form of f(z, y, y’, ---, y“””), the nature of the prob- 
lem, and the desired accuracy (all of which is a vast subject in itself). 


2. Comparison with Earlier Work. The chief novelty in the present arrangement 
is the systematic use of the u-series in terms of 5” to obtain (3) for any odd n (see 
also Milne below). Two other authors (Thiele, Buckingham), by employing the 
series for 1/u in terms of 8°”, give formulas for odd n that are closely related to (3), 
requiring just tabular arguments and avoiding the introduction of mid-interval 
arguments (as is done by Legendre, Oppolzer, Lindow). Presented chronologically, 
there is the following earlier work. 

Legendre [1] gives the symbolic formula for the (8/D)” series in 6°” and the 
first few coefficients up to n = 6. 

Oppolzer [2] gives the exact coefficients for (6/D) and (6/D)* up to 3”. His 
(8/D) coefficients checked with those in Salzer [5]. His (6/D)°* coefficients checked 
with A}, here, except for his coefficient Q.'* (= Ais) not in lowest terms by a 
factor of 9. 

Thiele [3] gives the exact values of the first five coefficients for D"" and (1/u)D~", 
which is the same as (5/D)” in terms of 6” and yo” up to m = 5, for n = 1(1)5. 

Lindow [4], who gives some central difference formulas up to triple quadrature, 
also gives the exact values of Aim , for m = 1(1)7. 

oA Salzer [5] tabulates the coefficients of 6/D, exactly through 8°”, then 18D through 
es 

Milne [6] happens to give 243m, m = 1(1)5, in the first of a series of formulas 

Zot+r 
ior [ , f(x)dx, r = 1(1)5, in terms of &'"fo . 
Zo-T 


Salzer [7] gives the coefficients of 5)°" and 6,°" obtained by k-fold quadrature of 
Everett’s formula; for k = 2, exactly up to m = 10, then 16D up to m = 24; for 
k = 3(1)6, exactly for m = 0 and 8S for m = 1(1)10. These differ from the other 
coefficients in that they occupy two lines for central differences instead of one. They 
are mentioned here because of their similar purpose and the large extent to which 
they have been tabulated. 

Buckingham [8] gives the coefficients of (6/D)" and (1/n)(6/D)", n = 1(1)4, 
through é°. As in Thiele [3], this includes an integration scheme involving just tabu- 
lar arguments for every n. Thus, by expressing (6/D)" for odd n as {(1/u)(6/D)"}u, 
and choosing 2» + h/2 for the argument, Buckingham obtains odd-order central 
differences of the integral, at mid-intervals, in terms of mean central even-order 
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differences, also at mid-intervals, so that both members involve y and f for just 
tabular arguments. However, it appears to the author that for n odd there is less 
total computation involved in using (3) for 4(é/D)”, where the slight extra work 
of finding ué"yo instead of 5”y:2 is more than compensated for by not having to 
average all the quantities sf, and &’"f, , as is done in the Buckingham-Thiele pro- 
cedure which uses (1/u)(4/D)” with the mean central differences yd°"f,2 . 


3. Integration Formulas for y” = f(x, y, y’,---, y 
ae + e B ss 4. 263 
180 1512 2 26800 149 68800 


eine 





n=1: py = (+e - 5” 


2 1 33787 2 1 57009 5! 
4 08648 24000 24 51889 44000 


_ 162 15071 6 4 26894 53969 538 
12504 63614 40000 99 78699 64291 20000 


2 68931 18531 *) 
ois fo 














4704 24411 73728 00000 
e a 31 289 317 
ee OS a, Oe 8 10 
oe Se RO+S 340 * 60480 36 28800° * 298 09000 ° 
_ 68 03477 32 + 32 03699 
261 53487 36000 627 68369 66400 
736 91749 is , 22 03877 95651 - 


~ 71137 48561 92000° * i0218 18843 43418 88000 


_ 15447 34732 56043 *) 
0 








g* 




















~ 337 20021 83332 82304 00000 

e 8s 4 661 
60480 57600 | 1596 67200 
_ 4 65967 2, 396079 gs 

52 30697 47200° ‘ 209 22789 88800 

95 95529 


16 
23712 49520 64000 ° + 2043 63768 68683 77600 


2143 27306 64071 *) f 
112 40007 27777 60768 00000 . 

> ee é 6 é 41 s, 491 10 
ditties ty =i (142 70 * 3021 7 25760° * 470 O0i00 
_3 41749) ow, SOOT 
17 43565 82400 13 07674 36800 
1704 03199 8 55137 58923 58 


2 13412 45685 76000 ° 5109 09421 71709 44000 


a 





n= 3: y= w+ + 5+ 





1 78574 25881 18 
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1721 38184 48999 i”) h 
~ 48 17145 97618 97472 00000 
a a 11 a 
- 5: 5 “ r® (1 7. a os S io. 
sginipe tte +3 + 33 9048 * 7 25760 ~ 708336 


BX 13283 3” + 5827 _ 
17 43565 82400 104 61394 94400 
_ 9 66067 gi 4. 4757 70541 - 
23712 49520 64000 364 93530 12264 96000 
_ 24 19396 16497 =) f 
6 88163 71088 42496 00000 7 























A ‘ ’, & Fy Fy 31 » 
am Sees 045+ 55+ 355 - stot eum! 
mi 27257 2 11581 _ 4 
3 11351 04000 6 22702 08000 
15 54079 is , 1 25853 54591 


= 18 
3908 65305 60000 si 1459 74120 49059 84000 , 


_ 150 48397 12643 ) j 
8 02857 66269 82912 00000 P 


General Dynamics/Astronautics 
San Diego, California 
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New Mersenne Primes 


By Alexander Hurwitz 


If p is prime, M, = 2” — 1 is called a Mersenne number. The primes M »; and 
Mz; were discovered by coding the Lucas-Lehmer test for the IBM 7090. These 
two new primes are the largest prime numbers known; for other large primes see 
Robinson [4]. The computing was done at the UCLA Compu ‘ng Facility. This test 
is described by the following theorem (see Lehmer [1, p. 443—.]) . 

Turorem. If S, = 4.and S,4; = S,” — 2 then M, is prime if and only if 8, = 0 
(mod M,). 

The test takes about 50 minutes of machine time for p = 4423. These results 
bring the number of known Mersenne primes to 20. The values of p for these twenty 
primes are listed in Table 1. 

If M, is prime it is of interest to know the sign of the least absolute penultimate 
residue, that is, whether S,.. = +2’ (mod M,) or S,. = —2’ (mod M,) where 
2r = p + 1. The Lucas-Lehmer test can also be used with S, = 10. The various 
penultimate residues of the known Mersenne primes were computed and the re- 
sults appear in Table 1 (see Robinson [3)). 

In addition to testing the above Mersenne primes each Mersenne number with 
p < 5000 was tested unless a factor of M, was known. The residues of S,_, (mod 
M,) are available for checking purposes. The results for 3300 < p < 5000 are sum- 
marized in Table 2. The computer program also found (see [3, p. 844]) that Mays) 
is not prime. 

The residue S,, (mod M,) for p > 3300 is output from the computer in a 
modified octal notation. That is, the residue is stored in the computer in 35 bit 
binary words and the output is a word by word conversion of the 35 bit words into 
octal (base 8) numbers. Thus the leading digit of each is quaternary (base 4). 
For p < 3300 the residue was printed in hexadecimal notation (see Robinson [3] 
and Riesel [2]). 











TABLE 1 
p j Si = 4 Si = 10 Pp Si = 4 Si = 10 
” Su. | a a ne 
2 | 107 - + 
3 + _ | 127 + 
5 + —- | 521 -— - 
7 _ —_ | 607 ~_ _ 
13 + 4 } 1279 _ _ 
17 ~ + | 2203 - ~ 
19 - + | 2281 = + 
31 + + 3217 - 
61 4 oe | 4253 + v 
89 _ + hel 


4423 -_ 


Received November 3, 1961. The preparation of this paper was sponsored by the U. S. 
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TABLE 2 











p R b R 
3301 72013 4241 11012 
3307 62061 4253 00000 
3313 10050 4259 46007 
3331 51270 4261 55632 
3343 76415 4283 74774 
3371 57040 4339 41356 
3373 36120 4349 74465 
3389 64705 4357 74271 
3413 50261 4363 61114 
3461 03241 4397 40174 
3463 57665 4409 51070 
3467 23046 4421 25131 
3469 21765 4423 00000 
3547 75574 4481 70216 
3559 45350 4493 36053 
3583 42507 4519 01571 
3607 45062 4523 22235 
3617 35431 4567 74267 
3631 14530 4583 46556 
3637 67413 4591 47243 
3643 04606 4621 74601 
3671 04031 4643 51444 
3673 01626 4651 00707 
3691 © 54715 4663 52442 
3697 53743 4673 40333 
3709 06427 4679 14305 
3739 22413 4703 54013 
3769 00747 4721 04420 
3821 52075 4729 40137 
3833 45453 4733 12774 
3847 57652 4783 77350 
3877 46507 4789 02364 
3881 34503 4799 04305 
3889 30737 4817 70020 
3919 16520 4831 33213 
3943 33442 4877 75412 
4007 17770 4889 24410 
4027 60265 4909 61113 
4049 31260 4937 26525 
4051 63236 4951 22971 
4091 55650 4973 03354 
4093 26670 4987 72275 
4111 20437 —_————— 
4133 66046 8191 03624 
4157 43640 

4159 62544 

4177 16076 

4201 53211 

4219 51756 

4231 51457 
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The five least significant octal digits of the residue appear in Table 2 for each 
p > 3300 tested. If p (3300 < p < 5000) is omitted from Table 2 a factor of 2” — 1 
is known. Some of these factors are not yet published but were communicated to 
the author by John Brillhart. 


My thanks to the Computing Facility for their help in this work, especially J. 
L. Selfridge and F. H. Hollander. 
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REVIEWS AND DESCRIPTIONS OF TABLES AND BOOKS 


18 [F].—Rocer Osporn, Tables of All Primitive Roots of Odd Primes Less than 
1000, University of Texas Press, Austin, 1961, 70 p., 30 cm. Price $3.00. 


This slim volume lists all 28,597 primitive roots of the 167 odd primes less than 
1000. These tables were computed on an IBM 650. The program and running times 
are not indicated. The most extensive earlier table, as noted by the author, is due 
to Chebyshev and extends to p = 353. 

There also is a small table of statistical information. Perhaps the most interest- 
ing column here lists the number of (positive) primitive roots less than p/2 for 
each prime p. Of the 87 primes = —1 (mod 4), eight have exactly one-half of 
their primitive roots less than p/2. The seven primes 223, 379, 462, 631, 691, 883, 
and 907 have more than one-half less than p/2. The remaining 72 primes have less 
than one-half there. The author associates this preponderance with the well-known 
fact that more than one-half of the quadratic residues of such primes lie in this 
interval. 

For the primes = +1 (mod 4) this column is clearly redundant, since it is 
easily seen that if g is a primitive root for such a prime then so is p — g. For these 
primes the real interval of interest is p/4 < g < 3p/4. Since the quadratic non- 
residues are in excess here, one would expect the primitive roots to also be pre- 
ponderantly in excess, since approximately three-fourths of all non-residues are 
primitive roots. 


D.S. 


19 [I, X].—D. 8S. Mrrrinovié & R. 8S. Mrrrinovié, Sur les nombres de Stirling et 
les nombres de Bernoulli de Vordre supérieur, Publ. Fac Elect. Univ. Belgrade 
(Série: Math. et Phys.), No. 43, 1960, 64 p. (French with Serbian summary.) 
The tables in this paper extend those given in previous papers, especially the 

three reviewed in Mathematics of Computation, v. 15, 1961, p. 107. The notation 

used is explained in that review. 

Table I (p. 15-44) gives (—)”"Cn* for k = 0(1)32, m = 33(1)50, and for k = 
33(1)49, m = k + 1(1)50, 

Table II (p. 45-50) gives SS,” ” for m = 33(1)49, n = m + 1(1)50, and also 
form = 50,2 = 51. 

Table III (p. 51-63) gives S,”” for m = 1(1)3, n = 201(1) 1000. 

The tables were computed on desk machines. Checks made by the authors were 
supplemented by comparison with Miksa’s unpublished tables and by many-figure 
computations made in laboratories at Liverpool, Rome, and Munich. A bibliog- 
raphy of 26 items is given. 

A. F. 


20 [K].—B. M. Bennett & P. Hsu, Significance Tests in a 2 X 2 Contingency 
Table: Further Extension of Finney-Latscha Tables, February 1961. Deposited 
in UMT File. 

These manuscript tables constitute an extension for A = 21(1)30 of tables pre- 
pared by Latscha for A = 16(1)20, and supersede the previous tables by the present 
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authors for A = 21(1)25. (See Review 9, Math. Comp., v. 15, 1961, p. 88-89.) The 
format and precision of those tables (four decimal places) is retained in this adden- 
dum. 


J. W. W. 


21 [K].—Courn R. Biytx & Davin W. Hutcuinson, Tables of Neyman Shortest 
Unbiased Confidence Intervals (a) for the Binomial Parameter (b) for the Poisson 
Parameter, (reproduced from Biometrika, v. 47, p. 381-391, v. 48, p. 191-194, 
respectively) University Press, London, 1960, 16 p., 28 em. Price 2s. 6d. 


Anscombe [1] observed that exact confidence intervals for a parameter in the 
distribution function of a discrete random variable could be obtained by adding 
to the sample value, X, of the discrete variable a randomly drawn value, Y, from 
the rectangular distribution on (0, 1). Eudey [2] has applied this idea in the case 
of the binomial parameter, p, to find the Neyman shortest unbiased confidence set. 
The present authors use Eudey’s equations for a uniformly most powerful level 
l-a test of p = p* vs p ~ p* based on an X in a sample of n, which give the ac- 
ceptance interval a(p*) determined by a value of Y in the form m + yo S X + 
Y <= m + in which m and n, are integers and 0 < y < 1,0 S 7:1 = 1. These 
are solved for yo and 7; in terms of m» and m, and the given X, n, and a. Then trial 
values of m) and n, are used until the resulting yy and 7; are both on (0, 1). The 
computation was carried out on the University of Illinois Digital Computer Lab- 
oratory’s ILLIAC. The program used for arbitrary n, a prints out mo + yo, m + 
for any equally spaced set of p* values. From these the Neyman shortest unbiased 
a-confidence set for p, X + Y € a(p*) can be read off to 2D. The tables give such 
95% and 99% confidence intervals for p to 2D for n = 2(1)24(2)50 and X + Y = 
0(.1)5.5 for n <= 10, 0(.1)1(.2)10 for 11 S nm S 19, 0(.1)1(.2)6(.5)15(1)17 for 
20 < n S 32, and 0(.2)2(.5)23(1)26 for 34 S n S 50. Forn, X + Y not tabled, 
one enters the table at n,n + 1 — (X + Y) and takes the reflection about p = 
3 of the interval given. 

Similar confidence intervals for the Poisson parameter, A, were found by the 
same method. The table gives Neyman shortest unbiased 95% confidence intervals 
for \ to 1D for X + Y = .01(.01).1(.02).2(.05) 1(.1) 10(.2)40(.5)55(1)59 and to 
the nearest integer for X + Y = 60(1)250. For the same values of X + Y, 99% 
confidence intervals are given to 1D for X + Y <S 54 and to the nearest integer 
for X + Y > 54. 

C. C. Craic 
The University of Michigan 
Ann Arbor, Michigan 

1. F. J. Anscomse, ‘The validity of comparative experiments,’’ J. Roy Statist. Soc. Ser. 

A. v. 111, 1948, p. 181- ‘211. 


2.M. W. Evpey, On the Treatment of a Discontinuous Random Variable, Technical 
Report No. 13 (1949), Statistical Laboratory, University of California, Berkeley. 


22 [L].—M. I. Zuurina & L. N. Karamazina, Tablifsy funkisit Lezhandra P_1)2+ir (x), 
Tom I (Tables of the Legendre functions P_;/24;, (x), Vol. 1), Izdatel’stov Akad. 
Nauk SSSR, Moscow, 1960, 320 p., 27 em., 2700 copies. Price 34.50 (now 37.95) 
rubles. 


This important volume belongs to the well-known series of Mathematical Tables 
of the Academy of Sciences of the USSR, and the tables were computed on the 
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high-speed electronic calculator STRELA at the Computational Center of the 
Academy. 

The Russian work has been concerned with the functions P_,2,:- (x), where r is 
real and z > —1. The functions are real, and satisfy the differential equation 


(1 — 2°)u” — 2eu’ — (4+ 7 )u = 0. 


The functions occur in potential problems relating, for example, to cones and hyper- 
boloids of revolution; they also occur in the Mehler-Fock inversion formulas [1]. 
The tables for —1 < x < 1 and z > 1 are given in Volumes I and II, respectively. 
The formulas given in the Introduction to Vol. I are limited to those which have 
some application in the range —1 < x < 1. The values were computed from 


P_spr4ie(x) ad F(3 ane tr, 3 + tr; 1; 4 ~* $2), 


where F'(a,b;c;z) denotes the hypergeometric function, and were checked by differ- 
encing. The main table (pages 13-312) gives values of P_s24:, (x) to 7S for r = 
0(0.01)50, x = +0.9(—0.1) —0.9, without differences. (It is stated that Vol. II, 
which the reviewer has not seen, gives values for x = 1.1(0.1)2(0.2)5(0.5) 10(10) 
60.) The interval in 7 has been made narrow because applications in mathematical 
physics frequently require integration with respect to 7. It'is stated that interpola- 
tion in rt may be performed by the three-point Lagrange formula with an error not 
exceeding 1.6 final units; it may be added that such an error can occur in only a 
small part of the table. Interpolation in z is naturally more troublesome, even well 
away from a logarithmic singularity ats = —1. 

An auxiliary table on pages 315-318 facilitates use of an asymptotic series for 
large 7; arc cos x and four coefficients which are functions of x are tabulated to 7D 
for x = 0.99( —0.01) —0.90, without differences. Values of the Bessel functions J 
and J, are required to be available for use with the auxiliary table. 

A useful bibliography of 16 items averages about one misprint per item in the 
five non-Russian titles, the most entertaining being MacRobert’s well-known book 
on “Spherical Harmonies” and a paper by Barnes on “Veneralized Legendre Func- 
tions.” 

The reviewer differenced about a hundred values without finding any error. 
Assuming its accuracy, this must be reckoned a valuable table. 


A. F. 


1. A. Erp&y1 et al, Higher Transcendental Functions, Vol. 1, McGraw-Hill, New York, 
1953, p. 175. 


23 [X].—-A. Cuarnes & W. W. Coorer, Management Models & Industrial Appli- 
cations of Linear Programming, v. 1, John Wiley & Sons, Inc., New York, 1961, 
xxiii + 471 p., 26 em. Price $11.95. 

This book is addressed to persons interested in the application of linear pro- 
gramming techniques to various aspects of management planning. Much of the 
material has been published previously by the authors in scattered journals and 
texts; however, this volume offers the advantage of a unified mathematical treat- 
ment of sundry topics in mathematical programming and managerial economics 
within the framework of adjacent-extreme-point techniques. 
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The earlier parts of this volume do not require mathematics beyond college 
algebra. The rudiments of linear programming theory and techniques are illustrated 
by means of simple numerical examples. An elementary machine loading problem 
is introduced to elucidate such concepts as linear model formulation, approximation 
of model types by scaling, the dual linear programming problem, and data accuracy 
and program sensitivity. The stepping-stone method for the classical Hitchcock 
transportation problem and transshipment problem are described at length. The 
procedure for dealing with degeneracy is also discussed. To explicate the concept 
of input-output analysis, a three-industry input-output model as an example of a 
“static, open Leontief model” is given. Feasible solutions are obtained by the Gauss 
elimination method. 

With the exception of the transportation algorithm, a rigorous mathematical 
treatment of the foregoing topics are presented in the succeeding parts of this vol- 
ume. Background material from the fields of matrix algebra, convex sets, and linear 
systems are developed and interpreted to provide an essentially self-contained 
account of the mathematics relevant to the managerial applications covered in the 
rest of the volume. 

Considerable attention is devoted to Dantzig’s simplex method for solving the 
general linear programming problem. The basic simplex algorithm is carefully 
explained and illustrated with the aid of numerical examples_and geometrical in- 
terpretations. Additional by-products and interpretations are obtained, such as the 
extension of the simplex calculations for analyzing the effects of altering (a) the 
stipulations vector, (b) the coefficients of the objective function, and (c) the strue- 
tural vectors. Also, the role of the simplex procedure as a tool for securing proofs 
of several important duality theorems in the field of linear inequalities is deftly 
portrayed. 

The application of delegation models to managerial economics is first examined 
along the lines of T. C. Koopman’s “activity analysis models.”” A major purpose of 
such models is the determination of rules which might be applied to guide the ac- 
tivities of a decentralized management organization. Koopman’s formulation is 
reduced to a series of special linear programming problems and their duals. “Effi- 
cient” solutions are obtained by the “spiral” method. Koopman’s concept of 
“efficiency” is then generalized to provide under certain circumstances more suitable 
criteria for managerial applications. 

Linear programming approaches to statistical problems involving inequality re- 
lationships are delineated and applied to a problem of determining an executive- 
compensation formula for an industrial concern. Moreover, the techniques employed 
to solve this problem provide an introduction to the use of adjacent-extreme-point 
methods to a variety of nonlinear problems encountered in management planning. 
Modifications of simplex criteria and procedures are developed for the case where 
a functional subject to linear constraints may be decomposed by linear transforma- 
tions into a sum of functionals involving only a single variable. The basic short- 
coming of this approach is that, in general, only a local optimum is guaranteed. 

A dynamic model for production scheduling at minimum cost when the costs 
are unknown is solved by means of “surrogate” techniques and “subhorizon” 
methods. Optimizing rules are enumerated and expounded for solving an actual 
example for which these methods were first devised. This is followed by a proof of 
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the optimizing properties of the rules. The effects of introducing costs, such as 
inventory charges, and additional constraints, such as storage limitations, are 
touched upon from the standpoint of possible variations in the length of sub- 
horizons. A generalized approach to this class of problem is explored via the Kuhn- 
Tucker theorem for nonlinear programming. 

The “classical” models of linear programming are presented with commendable 
clarity. Moreover, the adaptation of linear programming methods for solving non- 
linear types of management problems is aptly demonstrated. However, this re- 
viewer’s enthusiasm was tempered by the fact that the present edition abounds 
with errors resulting from an apparent cursory attempt at editing and proofreading. 
This reviewer recommends that the publishers prepare an errata sheet; otherwise, 
the intolerable number of typographical errors will vitiate the intrinsic merits of this 








book as a textbook and reference. 


Applied Mathematics Laboratory 


David Taylor Model Basin 
Washington 7, D. C. 


Mitton SIEGEL 


24 [X].—Roman Jaxosson, Editor, Proceedings of Symposia in Applied Mathe- 
matics, Vol. XII, “Structure of Language and its Mathematical Aspects,” 
American Mathematical Society, Providence, 1961, vi + 279 p., 26 cm. Price 


$7.80. 


Sponsored by the American Mathematical Society, the Association for Symbolic 
Logic, and the Linguistic Society of America, and cosponsored by the Institute for 
Defense Analyses under an Office of Naval Research contract, the symposium, held 
in April, 1960, included the following papers: 


W. V. Quine 
Noam Chomsky 
Hilary Putnam 
Henry Hiz 


Nelson Goodman 
Haskell B. Curry 
Yuen Ren Chao 


Murray Eden 
Morris Halle 


Robert Abernathy 
Hans G. Herzberger 
Anthony G. Oettinger 


Victor H. Yngve 

Gordon E. Peterson and 
Frank Harary 

Joachim Lambek 

H. A. Gleason, Jr. 

Benoit Mandelbrot 


Charles F. Hockett 
Rulon Wells 
Roman Jakobson 


Logic as a Source of Syntactical Insights 

On the Notion “Rule of Grammar” 

Some Issues in the Theory of Grammar 

Congrammaticality, Batteries of Transformations 
and Grammatical Categories 

Graphs for Linguistics 

Some Logical Aspects of Grammatical Structure 

Graphic and Phonetic Aspects of Linguistic and 
Mathematical Symbols 

On the Formalization of Handwriting 

On the Role of Simplicity in Linguistic Descrip- 
tions 

The Problem of Linguistic Equivalence 

The Joints of English 

Automatic Syntactic Analysis and the Pushdown 
Store 

The Depth Hypothesis 

Foundations in Phonemic Theory 


On the Calculus of Syntactic Types 

Genetic Relationship Among Languages 

On the Theory of Word Frequencies and on Re- 
lated Markovian Models of Discourse 

Grammar for the Hearer 

A Measure of Subjective Information 

Linguistics and Communication Theory 
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Some of the authors are concerned with preformal questions, i.e., with a discursive 
characterization of the substance of language; Quine, Putnam, Chao, Herzberger, 
and Jakobson seem to have such interests. Others are fully engaged with the con- 
struction of formal systems: Chomsky, Hiz, Curry, Halle, Abernathy, Peterson 
and Harary, Lambek, Mandelbrot, and Wells. Oettinger, Yngve, and Hockett aim 
at description of linguistic processors—natural or artificial—rather than at char- 
acterizations of language, although all three have formalisms to display. Eden, 
working on handwriting, might be placed with one of the latter two groups. Good- 
man’s contribution is the exposition of a branch of mathematics in its potential 
application to linguistic theory. Gleason shows the application of classification 
theory to a major branch of linguistics, the tracing of historical connections among 
languages. 

A cursory inspection of this volume would suggest that the “structure of lan- 
guage” is just its grammatical—or, more narrowly, syntactic—structure. Mandel- 
brot objects to the identification of “linguistics” and “grammar” (pp. 211-214), 
but mathematical formalization of linguistic theory is going forward more rapidly 
in syntax than in any other area, and it is, as Jakobson remarks (p. vi), mathematical 
logic and the theory of recursive functions in particular that is being applied. 
Mandelbrot seems to agree with his opponents that “statistical” and “grammatical” 
models are ‘“contradictory.”” He supposes that they must remain so; a different 
possibility is that grammatical models will furnish a structure on which statistical 
models can be developed. Grammar in any case is not the whole of linguistics, and 
problems like Gleason’s will probably be brought to computing centers more often 
in the future. 

Computational linguistics has been hampered by lack of sufficient and sufficiently 
sound publications in mathematical linguistics; this volume should be studied by 
any linguist or mathematician who proposes to program syntactic operations, 
whether for research purposes or in connection with such applications as machine 
translation. 

Davin G. Hays 
The RAND Corporation 
Santa Monica, California 


25 [Z].—Dona.p P. Eckman, Editor, Systems: Research & Design, John Wiley & 
Sons, Inc., New York, 1961, xiii + 310 p., 23 cm. Price $8.50. 


This book is the Proceedings of the First Systems Symposium at Case Institute 
of Technology. It contains a Foreword, a Preface, and fourteen papers concerning 
systems research and systems design. The fourteen papers vary in style, most 
noticeably with regard to bibliographic reference. Some are simply advice from the 
author without reference to other work, others have extensive bibliographies. Only 
one pertains directly to the mathematics of computation, “‘A problem in the design 
of large-scale digital computer systems” by R. J. Nelson. This paper is devoted 
almost entirely to the problem of designing a machine which would be efficient in 
selecting the largest number of a set and (by implication) in other sorting problems. 
No specific design is arrived at, but a facility for scanning a region of the memory 
is suggested; the ideas may mislead some readers if they are unfamiliar with thresh- 
old search commands such as that of the Control Data Corporation 1604 computer 
and with the engineering details of comparison circuits. 
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Other papers have implications connected with the mathematics of computation, 
as would be expected in any current book on large systems. Thus in the Foreword, 
Simon Ramo remarks that “it could be said that systems engineering in today’s 
sense became possible only with the introduction of the large digital computer.” 
However, the papers in this volume contribute few direct suggestions concerning 
this use, and concern themselves largely with other general and specific aspects of 
systems engineering. 


C. B. Tompkins 
University of California 


Los Angeles, California 


26 [Z).—Daniet D. McCracken, A Guide to FORTRAN Programming, John 
Wiley & Sons, Inc., New York, 1961, viii + 88 p., 28 em. Price $2.95. 


The usefulness of Fortran as an automatic programming system available on 
many different computers has prompted Dr. McCracken to publish this guide. It 
is addressed to people who have no programming experience but have a requirement 
to accomplish scientific computation or wish to get some appreciation of how this 
can be done. 

The guide is developed pedagogically, with numerous examples, and includes a 
set of detailed case studies which provide examples from several fields of effort. 
These case studies illustrate the essential features of Fortran and suggest the range 
of its applicability. 

An appendix summarizes the characteristics of a number of Fortran systems that 
have been established for different computers. 


ABRAHAM SINKOV 
National Security Agency 


Department of Defense 
Washington 25, D. C. 


27 [Z].—Francis J. Murray, Mathematical Machines, Vol. 1 and 2, Columbia 
University Press, New York, 1961. V. 1, vii + 300 p., 26 em. Price $12.50. V. 
2, vii + 365 p., 26 em. Price $17.50. 


Volume one of Professor Murray’s two-volume work on mathematical machines, 
is concerned with digital computers. There are two parts in Volume 1: part I on 
desk calculators and punched card machines, and part II on automatic sequence 
digital calculators. These digital devices are presented in the order of increasing 
competence and complexity. 

In part I, there are eight chapters. The first four chapters describe desk calculators, 
from the basic idea of register and counter to the description of many commercial 
automatic calculators. Chapter 5 covers electrical counters and accumulators. 
Punched card machines are presented in Chapters 6 and 7, and sequence calculators 
such as calculating punch and electronic calculator in Chapter 8. 

Part II consists of ten chapters. The first four chapters describe the logic aspect 
of the computer as well as digital arithmetic. Chapter 5 is a general discussion on 
the use of Boolean analysis. Chapter 6 is concerned with circuit elements. The pro- 
gramming aspects are covered in Chapters 7, 8, and 9. Chapter 10 is a very brief 
survey of digital computers. 
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In this volume, the author succeeded in many cases in bringing out the principles 
and fundamental ideas. An example is the exposition on desk calculators. Although 
the material is mostly descriptive, it will serve a useful purpose as a general refer- 
ence. 

Volume two of Professor Murray’s work on mathematical machines presents the 
subject of analog devices. There are three parts: part III on continuous computers, 
part IV on true analogs, and part V on mathematical instruments. 

Part III consists of fifteen chapters. After a brief introduction in Chapter 1, 
Professor Murray describes mechanical adders, multipliers, dividers, and other 
mechanical components in Chapter 2. Cams, gears, and their computing applica- 
tions constitute Chapter 3. This is followed by an excellent presentation on mechani- 
cal integrators, differentiators, and amplifiers in Chapter 4. Chapter 5 is a review of 
circuit theory. Computation by using potentiometers and condensers are described 
in Chapter 6, vacuum tube amplifiers in Chapter 7, electromechanical components 
of D’Arsonval movement, watt-hour meters, and synchros in Chapter 8, electrical 
multipliers including time division multipliers, strain gauge multipliers, step multi- 
pliers, cathode ray multipliers in Chapter 9, and function generation by using 
mechanical, electromechanical and electronic means in Chapter 10. Chapters 11 
through 13 describe equation solution: linear equations in Chapter 11, harmonic 
analysis and polynomial equations in Chapter 12, differential equations in Chapter 
13, and error analysis in Chapter 14. Chapter 15, the last chapter of this part, dis- 
cusses the use of digital check solutions obtained by using numerical methods when 
the analog solution has narrowed down the range of parameters. 

Part IV, consisting of nine chapters, presents the idea of true analogs. True ana- 
logs are direct analogies on which measurements can be taken more conveniently 
or more economically than the analog devices described in part III. The author 
examines the theory of true analogs and includes descriptions of dimension theory, 
models, and principles of spatial relationships. True analogs that are described 
include the use of electrolytic tanks, electrically conductive sheets, stretched mem- 
branes, photoelastic models, and electromechanical analogies. 

Part V consists of five chapters. It deals with mathematical instruments that 
operate on data in a specified form and perform a few mathematical operations. 
These devices include slide rule, plotting devices, planimeters, integrometers, 
integraphs, and other geometrical and trigonometrical devices. This part is rather 
unique. 

This volume again emphasizes principles. A significant portion describes me- 
chanical analog devices. The treatment of analog devices in volume two is more ex- 
tensive than that of digital computers. 

As mentioned in the book, this work was sponsored by the Office of Naval Re- 
search. These two volumes are a contribution to the study of mathematical ma- 
chines, and Columbia University Press deserves credit for an excellent printing 
job. 


Yaouan Cuvu 
Electrical Engineering Department 
University of Maryland 
College Park, Maryland 
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28 [Z].—E. L. Witiey, Marion Triss, A. D’Acapryerr, B. J. Gissens & MI- 
CHELLE CLARK, Some Commercial Autocodes, Academic Press, Inc., New York, 
1961, vii + 53 p., 23 em. Price $2.50. 


Some Commerical Autocodes is a study of nine programming languages applicable 
to commercial data processing problems, compiled in a tabular form by language 
elements. The study is based upon information available in December 1960 and 
does not represent the final specifications for some languages which have been, or 
are being, implemented for the various computers. 


Frances E. HoLBerRTON 
Applied Mathematics Laboratory 


David Taylor Model Basin 
Washington 7, D. C. 








TABLE ERRATA 


308.—A. Erpéiy1, W. Maenus, F. Ospernerrincer & F. Tricomi, Higher Tran- 
scendental Functions, McGraw-Hill Book Co., Inc., New York, 1953. 


The following corrections should be made in this work: 
Volume I 
P. 104, eq. (43); for (ec — a)F(c + 1) read (c — a)zF(e + 1). 
P. 145, eq. (24): replace italic P and Q by their roman equivalents. 
P. 150, second of eqs. (13): for i, read -2. 
Volume IT 
P. 321, eq. (22): for k’, read k”; and for E(0, k), read E(0, k’). 


J. C. Cooke 


Aerodynamics Department 
Royal Aircraft Establishment 
Farnborough, Hampshire 
England 


309.—MervIN E. Mutter, “An inverse method for the generation of random nor- 
mal deviates on large-scale computers,”’ MT AC, v. 12, 1958, p. 167-174. 


The following errors have been noted in Table 5, “Inverse Values for the Nor- 
mal Distribution”’: 


j F(x;) xj 
reads should read 

36 0.64062 500 0.36013 003 0.36012 989 

92 0.85937 500 1.07750 557 1.07751 557 

96 0.87500 000 1.15035 938 1.15034 938 
100 0.89062 500 1.22984 876 1.22985 876 
102 0.89843 750 1.27268 865 1.27269 865 
110 0.92968 750 1.47345 903 1.47346 759 
116 0.95312 5 1.67594 192 1.67593 973 
119 0.96484 375 1.80989 233 1.80989 224 


G. Miniter CLarK 
Ohio State University 
Columbus 12, Ohio 


310.—D. J. Finney, “The Fisher-Yates test of significance in 2 X 2 contingency 
tables,” Biometrika, v. 35, Parts 1 and 2, May 1948. 


These tables have been checked against Tables of the Hypergeometric Probabil- 
ity Distribution, by G. J. Lieberman and D. B. Owen, Stanford University Press, 
1961. All the entries were found to be correct, except for the following typographical 
error: 

p. 149 A=6,B=5,a=6 Probability = 0.025 
for 90 015 read 1 O15. 
261 
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This error is reproduced in Table 38 on page 188 of Biometrika Tables for Statis- 
ticians, Volume 1, by E. 8S. Pearson and H. O. Hartley, University Press, Cambridge, 
1954. 

Anna M. GLINsKI 


JoHn VAN Dyke 
National Bureau of Standards 
Washington 25, D. C. 


311.—R. Latscna, “Tests of significance in a 2 X 2 contingency table: extension 
of Finney’s table,” Biometrika, v. 40, Parts 1 and 2, June 1953, p. 74-86. 
These tables have been checked against the Lieberman-Owen Tables of the 

Hypergeometric Probability Distribution, and the following errors noted. 


A B a prob. for read 

16 10 14 0.05 4 .018 4 .017 
16 10 14 0.025 4 .018 4 .017 
16 4 15 0.005 1 .001 0 .001 
17 4 16 0.05 , On 1 .012 
17 4 16 0.025 ; oer 1 .O12 
19 16 13 0.025 4 .012 4 .012 
19 8 15 0.05 2 .013 2 .014 
19 8 15 0.025 2 .013 2 .014 
19 6 19 0.05 4 .050— 4 .050 
20 15 17 0.005 5 .002 5 .003 
20 12 19 0.05 7 .019 7 .018 
20 12 19 0.025 7 .O19 7 .018 


In order to be consistent with the method of construction for this table, in which 
the value of b recorded is the greatest significant value for which the corresponding 
probability is less than or equal to the probability shown at the head of the column, 
the following additional line should be inserted in the appropriate place in the 
table: 


Probability 
A B a 0.05 0.025 0.01 0.005 
19 1 19 0 .050 --- —---- -_—— 
Anna M. GLINSKI 
JOHN VAN DYKE 
Corrigenda 


AnpREs Zavrotsky, ‘“Construccion de una escala continua de las operaciones 
aritmeticas,” Math. Comp., Review 63, v. 15, 1961, p. 299-300. 


On page 300, line 7, instead of L” x = H(Gx — 1), read L" x = H(Gr — n). 
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R. T. Ostrowski & K. D. Van Doren, “Ona theorem of Mann on latin squares,” 
Math. Comp., v. 15, 1961, p. 293-295. 
PUTTS 1/10\" 1/10)’ 
On page 294, line 18 from the bottom, for 7 : 2 ag 15,876, read = “) Bw 15,876. 


/ 


ARNOLD N. Lowan, “On th numerical treatinent of heat conduction problems 
with mixed boundary conditions,” Math. Comp., v. 14, 1960, p. 266-270. 


For equations (13), (14), and (15) on page 269, read 
Th an+1 = aT 11,0 + (1 es 2a oo B) Fite + aT nassin + BT 2.0 + Una.n 


a/Arsh<M (13) 

T wend = BT a pan + OT yin + (1 — @ — 28)T aan + BT a asin (14) 
+ Uns. l1<k<WN 

Trwngt = BT awa + aTiawan + (1 — 2a — B)Tiwn (15) 
+ aPaase + Us.n0 eo/Ac Sh<M 


where U;,;,., and Uy.» and U;.y,, are the same as previously given. In addition, 
for points bounded on two sides by heat fluxes, the equations must be further modi- 
fied to give 


Tuan = @Tyian + (1 — @ — B)T wan + BT usin + Unan 
+ Uyin for h= M, k= 1 

and 

Tune = BT uwan + aT uawn + (lL — @ — B)T awn + Unis 


+Usv. for h=M, k=N 


Norman J. McCormick 
Ann Arbor, Michigan 
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